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ABSTRACT 


Propulsion  losses  are  increased  by  added  drag  due  to 
steering  of  the  ship.  A  carefully  designed  automatic 
steering  control  provides  the  desired  heading  while  it 
simultaneously  minimizes  the  rudder  activity  and  holds  the 
potential  for  reducing  propulsive  losses. 

A  computer  model  of  the  SL-7  containership  along  with  a 
cascaded  controller  (one  pole,  one  zero)  were  coupled  to  a 
function  minimization  subroutine  and  a  sea  state  generator 
program.  This  scheme  provided  the  appropriate  controller 
parameters  in  order  to  accomplish  the  best  performance. 

The  model  was  tested  in  calm  waters  and  sea  states 
(regular  and  irregular)  as  well,  for  a  certain  speed  and 
different  encounter  wave  angles  and  encounter  frequencies. 

Also,  an  adaptive  control  was  studied  which  updates  the 
controller  parameters  while  either  the  environmental  condi¬ 
tions  or  the  ship's  steering  characteristics  change  in  order 
to  maintain  optimal  steering  performance. 
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I .  INTRODUCTION 


The  economics  associated  with  ship  operations  have 
necessitated  an  examination  of  the  losses  associated  with 
the  motion  of  an  automatically  steered  ship  in  a  seaway. 

Four  major  areas  where  fuel  losses  occur  during  the 
operation  of  a  ship  have  been  identified  [Ref.  1,  2].  These 
areas  on  existing  s t earn / dies e 1  tankers  are  shown  below: 

•  Power  plant  and  auxiliaries 

•  Propeller  efficiency 

•  Hull  resistance 

•  Steering  and  navigation 

An  optimized  autopilot  design  would  provide  effective 
steering  control  with  associated  cost  savings  due  to 
reducing  fuel  consumption. 

An  appropriate  computer  model  which  represents  the  ship 
is  necessary  for  studies  leading  to  appropriate  controller 
design.  Chapter  2  introduces  the  development  of  two  of  these 
models . 

Chapter  3  addresses  the  formulation  of  a  performance 
criterion  which  represents  the  added  drag  due  to  steering  of 
the  ship. 

Using  the  equations  of  motion  as  a  model  of  the  ship  and 
a  function  minimization  subroutine  we  proceed  to  the 
controller  design  for  regular  seas  (deterministic  model  for 
the  seaway)  in  Chapter  4,  and  for  irregular  seas  (nondeter- 
ministic  model)  in  Chapter  5.  The  function  minimization 
subroutine  used  was  BOXPLX  and  was  programmed  by  R.  R. 
Hilleary  of  the  Naval  Postgraduate  School  Computer  Center 
[Ref.  3].  It  will  find  the  minimum  of  any  arbitrary  func¬ 
tion,  linear  or  nonlinear,  subject  to  explicit  constraints 
of  the  variables  or  implicit  constraints  on  functions  of  the 
variables . 

11 


Chapter  6  introduces  another  function  minimisation 
subroutine  appropriate  for  onboard  use. 

An  adaptive  control,  which  updates  the  controller  param¬ 
eters  when  the  environmental  conditions  or  the  ship's  course 
change,  is  studied  in  Chapter  7. 

Conclusions  drawn  from  these  experiments  and  recommenda¬ 
tions  for  future  studies  are  addressed  in  Chapter  8. 
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II.  COMPUTER  MODELS  OF  THE  SHIP 


A  nontrivial  part  of  any  control  problem  is  modelling 
the  process.  Thus,  an  appropriate  computer  model  which 
represents  the  ship  is  necessary.  The  best  representation  of 
the  ship’s  steering  dynamics  is  a  Taylor's  series  expansion 
of  the  force  and  moment  relationships  around  a  selected 
steady  state  operating  point.  The  equations  obtained  in  this 
way  are  known  as  the  equations  of  motion  [Ref.  4],  and  the 
formulation  in  the  computer  program  is  indicated  in  Appendix 
A.  This  computer  program  was  developed  by  using  known  avail¬ 
able  data  for  the  SL-7  containership  and  by  implementation 
of  the  scheme  in  Figure  2.1  [Ref.  5]. 

In  this  scheme  the  function  minimisation  subroutine  is 
fed  by  the  yaw  error  1 1/  and  rudder  angle  (5  ,  computes  the 
performance  criterion  J  and  adjusts  the  controller  free 
parameters  in  order  to  minimize  J. 

A  second  model  for  the  ship- steering  dynamics  represen¬ 
tation  is  the  Nomoto  model.  Figure  2.2  indicates  the  second 
and  third  order  Nomoto  transfer  functions  while  Figure  2.3 
indicates  the  appropriate  scheme  used  for  obtaining  these 
models  from  the  equations  of  motion.  Appendix  A  includes 
the  computer  program  used  for  the  Nomoto  third  order  model 
determination . 


A  yaw  command  is  applied  as  input  in  the  scheme  in 
Figure  2.3  and  the  difference  of  the  signals  l b  and  t//  is 

r  a  .  -  ...  M  EQ 

fed  to  the  function  minimization  subroutine  which  attempts 
to  adjust  the  free  parameters  of  the  Nomoto  plant  in  order 
to  minimize  the  performance  criterion  J. 

Simulation  runs  indicate  that  the  resulting  Nomoto 
models  are  obtained  with  resulting  J  close  to  zero. 
However  ,  in  this  study  the  equations  of  motion 


Figure  2.1  Controller  Optimization  Scheme 
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Figure  2.2  Nomoto  Transfer  Functions 
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Figure  2.3  Nomoto  Model  Determination  Scheme 
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representation  was  adopted  because  the  system  is  dynamic 
and  use  of  the  N'oinoto  model  representation  implies  addi¬ 
tional  computer  use.  On  the  other  hand,  frequency  domain 
studies  were  carried  out  using  the  Nomoto  representation 
since  this  represent  it  ion  is  easier  handled. 


AN  ADEQUATE  PERFORMANCE  CRITERION 


A.  CRITERION  BASED  ON  TRUE  ADDED  RESISTANCE 

The  performance  criterion  which  characterizes  propulsion 
losses  due  to  steering  may  be  shown  to  be  that  derived  from 
excess  power  consumption  per  unit  distance  caused  due  to 
steering  [Ref.  1,  6].  The  added  resistance  due  to  steering 
can  be  related  to  the  surge  or  thrust  equation  where  the 
total  instantaneous  surge  relevant  to  steering  is 

AX=[m>(  p/2)LAX^]yr*l/2[(p  /2)AX;y]v2  (3.1) 

♦L/2[(P  /2)AX^  Uj]  5  J 


where 


m  =  mass  of  ship 
P  =  density  of  sea  water 

L  =  ship's  length  between  perpendiculars 
A  =  L* 

U  =  ship's  water  speed 
v  =  sway  velocity 
r  =  yaw  rate  of  ship 
<5  =  rudder  angle 

Xyp  =  force  coefficient  due  to  yaw/sway  (posi¬ 
tive  ) 

Xk  =  f  orce  coefficient  due  to  rudder  angle 
(negative ) 

X’  =  force  coefficient  due  to  sway 


Since  the  sway  velocity  of  the  ship  is  small  we  can 
neglect  the  term  which  includes  the  square  of  the  sway 
velocity  in  the  previous  equation.  From  this  the  mean  surge 
relevant  to  steering  may  be  written  as 


AX=[m*(  p/2)LAX^r]  (u0rfl  2)cos(  ^  -  <p  ) 
*  [  (  P  /  2  )  AX  U 2  ]  (  5 2  /  2  )  *  r 


(3.2) 


Where  =  amplitude  of  sway  velocity 

=  amplitude  of  yaw  rate 

(5 0  =  amplitude  of  rudder  angle 

( P~{D=  phase  difference  between  sway  and 
V  r 

yaw  rate 


A  performance  criterion  for  added  resistance  due  to 
steering  may  be  formulated  as 


J=  lim  (l/2T)/(  -  vr  +  'Y  U2  (5  2  )dt 

t~*o°  Jn 


(3.3) 


where  (7  and  ~y  are  constants 


Accurate  knowledge  of  the  nonlinear  coefficients  X’  ’  and 
Xjj  is  required  for  the  accuracy  of  such  a  criterion.  In 
addition  the  criterion  itself  suffers  from  the  disadvantage 
that  sway  velocity  measurements  are  not  available. 
Normalizing  the  last  equation  the  performance  criterion  will 


W^1/2T)jf(_Xvr+52)dt 


(3.4) 


where  A= { 2 [m+ ( p  / 2 )LA] X  ’  }  /  [  ( p /2)X£S  U2] 


Table  I  indicates  the  values  of  A.  for  the  operating 


range  of  speed  of  the  ship  studied. 
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TABLE  I 

Weighting  factor 


Ship's  speed  (knots) 

23 

32 


21.5330 
10.4215 
5 .3900 


B.  CRITERION  BASED  ON  APPROXIMATE  ADDED  RESISTANCE 

Empirical  criteria  based  on  an  approximation  to  added 
resistance  may  also  be  derived.  A  semiempirical  criterion 
for  measuring  the  relative  performance  was  developed 
[Ref.  7],  based  on  the  assumption  of  small  amplitude  oscil¬ 
lations  around  the  steady-state  pivot  point  of  the  ship 
during  yawing  at  the  ship/steering  system  natural  frequency. 
This  may  be  extended  and  an  alternative  criterion  for  added 
resistance  will  be 


J= lim  (1/2T 

7—  oo 


)J cAi +  52)dt 


(3.5) 


where 


X  =X(x)=  [2m(l*Xyr)(0P/LXJ]/[(  p/2)LXg(j-U2] 
C(P /2)LAX;r]/m 

5P  =  distance  from  center  of  gravity  to  pivot 


center 


OJ  =  natural  frequency  (closed  loop  ship 
steering  control) 

P  -  ship's  perturbation  yaw  angle 


The  values  of  A  as  a  function  of  ship's  speed  are  given 
by  Table  II. 

A  closed  loop  system  natural  frequency  (x)  of  around  0.05 
rads  per  sec  has  the  potential  to  attenuate  the  effects  of 


TABLE  II 

Weighting  factor 


Ship  ’  s 


s 

15 
23 
32 


eed 


(knots  ) 


X 


x 


720 

251 

680 


seaway  disturbance  in  the  range  of  encounter  angles  where 
added  resistance  due  to  steering  is  important  [Ref.  6].  The 
weighting  factor  for  the  operating  range  of  the  ship  is 
shown  in  Table  III. 


TABLE  III 
Weighting  factor 


Ship '  s 


s^eed 

23 

32 


(knots ) 


X 

X 

16.796 


Equation  3.5  is  used  as  a  performance 
study.  It  is  an  approximation  but  it 
onboard  use  since  ship’s  perturbation 
rudder  angle  5  are  measurable. 

C.  WEIGHTING  FACTOR  STUDY. 

The  weighting  factor  X  given  by  Table  III  used  in  equa¬ 
tion  3.5,  plays  an  important  role  in  terms  of  the  optimal 
controller  parameters  determination.  Some  investigation  is 
necessary  in  order  to  verify  the  accuracy  of  the  results, 
since  the  values  of  X  of  Table  III  are  determined  based  on 
the  assumption  that  the  closed  loop  system's  natural 
frequency  is  around  0.05  rads  per  sec  [Ref.  1].  Frequency 
domain  techniques  were  used  for  this  purpose.  Using  the 
Nomoto  third  order  model  representation  of  the  ship  and 
available  controller  parameters  from  Chapter  4  for  sea  state 


criterion  for  this 
is  convenient  for 
yaw  angle  \L  and 


4,  encounter  frequency  1.5  rads  per  sec,  encounter  angle 
150°  and  ship’s  spc  ad  23  knots  we  found  that  the  closed  loop 
bandwidth  of  the  system  is  0.04  rads  per  sec,  as  indicated 
in  Figure  3.1,  which  is  not  close  enough  to  0.05  rads  per 
sec . 

For  the  same  sea  conditions  and  ship  speed,  with  the 
assumption  that  the  closed  loop  natural  frequency  of  the 
system  is  not  0.05  rads  per  sec  but  0.04  rads  per  sec,  a  new 
value  X=5.734  was  obtained  and  the  frequency  domain  tech¬ 
niques  result  in  a  new  bandwidth  0.035  rads  per  sec  as  is 
indicated  in  Figure  3.2. 

Clearly,  the  values  of  X  given  by  Table  III  and  used  in 
this  study  are  not  the  best.  Unfortunately,  since  the  full 
hydrodynamic  coefficients  of  the  SL-7  containership  are  not 
known  we  can’t  develop  the  surge  equation  and  thus  it  is 
still  impossible  to  determine  accurate  values  for  the 
weighting  factor  X. 
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IV.  REGULAR  SEAS 


CONTROLLER  DESIGN 


We  have  already  defined  a  suitable  and  sufficiently 
accurate  ship  computer  model  and  the  system's  performance 
criterion,  as  well.  The  remaining  task  is  to  determine  a 
representation  of  the  external  disturbances  imparted  to  the 
ship  by  the  sea,  before  the  system's  performance  in  a  seaway 
can  be  evaluated.  A  correct  model  of  the  seaway  itself  is 
essential  to  representative  modeling  of  forces  and  moments 
exerted  on  the  ship  by  it. 

At  this  point  we  will  use  the  regular  sea  model  as  sea 
representation.  The  properties  of  regular  seas  are  well 
defined.  The  wave  crests  are  assumed  to  be  straight,  infi¬ 
nitely  long,  parallel  and  equally  spaced  with  constant  wave 
height.  The  waves  progress  in  a  direction  perpendicular  to 
the  crest  line  at  a  uniform  velocity.  However  the  sea  is 
never  regular.  It  is  a  random  phenomenon  where  waves  are 
continually  changing  in  height,  length  and  breadth  [Ref.  8]. 

The  forces  exerted  by  the  regular  sea  have  the  form 


F=6J(7R/-cos(  (Jet*  $y  ) 


(4.1) 


where  C0^=  significant  wave  height 


Ry =  exciting  force 
( jJq =  encounter  frequency 
■$•=  phase  angle 

The  correspondence  between  sea  state  and  wave  height  is 
indicated  in  Table  IV  [Ref.  9], 

The  exciting  forces  Ry  for  different  encounter  frequen¬ 
cies  and  encounter  angles  were  obtained  from  the  sea  state 
generator  program  [Ref.  10]. 
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TABLE  IV 

Sea  state  vs  Wave  height 


Sea  state 
(Beaufort  scale) 
0 
1 
2 

3 

4 

5 

6 

7 

8 
9 


Range  for  wave  height 
(Feet) 

0 . 0 
0.32 

0.65-0.98 

1.96-3.28 

3.28-4.92 

6.56-8.20 

9.84-13.1 

13. 1- 18.2 

18.2- 24.6 
23.0-32.9 


Appendix  B  indicates  the  regular  seastate  formulation  in 
the  FORTRAN  program  used  for  obtaining  the  controller  param¬ 
eters  . 

The  controller  used  in  the  entire  study  has  one  pole-one 
zero  and  the  form  is  indicated  in  Figure  4.1.  This 
controller  seems  to  have  the  best  performance  in  calm  waters 
and  in  seaway  [Ref.  5], 


Figure  4.1  Controller  Used  in  this  Study 

The  optimized  controller  parameters  and  the  cost  J  for 
23  knots  speed,  sea  states  4-6-7-9,  different  encounter 
angles  and  various  encounter  frequencies  are  indicated  in 
Tables  V, VI, VII  and  VIII. 
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Studying  the  Tables  V  through  VIII  we  can  draw  the 
following  conclusions: 

•  For  a  particular  encounter  angle  and  encounter 
frequency  the  higher  the  sea  state  the  higher  the  cost. 

•  For  the  same  sea  state  the  cost  becomes  smaller  for 
higher  encounter  frequencies. 

•  For  encounter  frequency  0.2  rads  per  sec  the  maximum 
cost  occurs  at  60°  encounter  angle  for  all  tested  sea 
states . 

•  For  0.6  and  0.75  rads  per  sec  encounter  frequency  the 
maximum  cost  occurs  at  120°  encounter  angle  for  all 
tested  sea  states. 

•  For  1.5  rads  per  sec  encounter  frequency  the  maximum 
cost  occurs  at  90°  encounter  angle  for  all  tested  sea 
states . 

•  For  0.4  rads  per  sec  encounter  frequency  the  maximum 
cost  occurs  at  90°  encounter  angle  for  sea  states  4,6 
and  7  while  at  sea  state  9  the  maximum  cost  occurs  at 
60a  encounter  angle. 

Appendix  C  provides  the  computer  program  necessary  to 
achieve  the  system's  response.  Some  typical  responses  are 
indicated  in  Figures  4.2,  4.3,  4.4,  4.5,  4.6,  4.7,  4.8,  4.9. 
It  is  obvious  that  as  the  sea  state  goes  to  heavier  seas  the 
rudder  and  yaw  perturbations  become  larger. 

An  attempt  to  determine  how  accurate  the  controller 
parameters  must  be  for  a  particular  situation,  leads  to  the 
conclusion  that  high  accuracy  isn't  required.  Keeping  two 
parameters  fixed  each  time  and  vary  the  third  we  can  see 
(Figures  4.10,  4.11,  4.12,  4.13,  4.14)  that  the  cost  doesn't 
change  appreciably  in  the  vicinity  of  the  actual  value. 

Figures  4.2,  4.3,  4.4,  4.5,  4.6,  4.7,  4.8,  4.9  indicate 

that  the  yaw  and  rudder  excursions  are  less  than  1°.  This 
just  seems  strange,  though  it  may  be  because  of  the  opti¬ 
mization  of  the  filter.  We  tried  to  investigate  that  by 
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using  the  optimal  filter  for  sea  state  9,  encounter 
frequency  1.5  rads  per  sec,  encounter  angle  120°  and  run  it 
in  sea  state  4,  keeping  the  same  encounter  frequency  and 
angle.  The  yaw  and  rudder  excursions,  even  if  they  became 
larger,  remained  less  than  la .  The  parameters  of  those  two 
filters  are  close  and  the  reason  might  be  the  flatness  of 
the  cost  surface.  Second  attempt  led  to  more  interesting 
results.  Using  the  same  filter  and  run  it  in  sea  state  4, 
encounter  frequency  0.4  rads  per  sec  and  encounter  angle 
060°,  the  system  becomes  unstable  (Figures  4.15,  4.16). 


TABLE  V 


Optimal  Controller  Parameters 

Sea  State  4 


for  Regular  Sea 


Encounter  Frequency  0.2  rads 


encounter 
angle (degrees ) 


K 1 


T1 


per  sec 
T2 


0 

0.5221067 

66.3312231 

12.8332741 

30 

1.0488701 

61.9309387 

15 . 9266357 

60 

1.2036362 

54.5'533295 

16.0245972 

90 

1.3178062 

49 . 7426453 

14.8329315 

120 

1.3984699 

46 .9797058 

13 . 9525757 

150 

1.4502153 

45.4263306 

13 . 3351599 

180 

0.7195223 

25.2119598 

14. 1219782 

33 


609E- 
4819 
612223 
703524 
355578 
3567196 
.  345E-28 


Encounter 

encounter  Kl 

angle (degrees ) 


F requencv  0 


4  rads 


per 

T2 


sec 


0 

0.5221067 

66  .  3312231 

12 .8332741 

0.517E-35 

30 

0.7234516 

74.7846533 

47 . 9879893 

0.54673 

60 

0.8730211 

92.8420868 

50.0454407 

1 .  194746 

90 

2.8910360 

28.9679871 

10.0176315 

2.536161 

120 

2.6232796 

27.9702454 

8 .8327761 

0.2349457 

150 

2.6408129 

27.7937927 

8.6838455 

0.0545772 

180 

0.7195223 

25.2119598 

14.1219782 

0.536E-32 

Encounter  F 

requencv  0 . 6 

rads  per  sec 

eUCUUUL e  L 

angle (degrees ) 


0 

0.5221067 

.30 

2.0506487 

60 

1.0504951 

90 

2.1496201 

120 

2.1221027 

150 

1.9617786 

180 

0.7195223 

Encounter  F 

66.3312231 

1.2107763 

0.2329917 

1.2607498 

0.9715905 

0.8480816 

25.2119598 


12 

5 

19 

5 

6 
8 

14 


8332741 

5998001 

0010876 

3318434 

9064713 

3953667 

1219782 


0 

0 

0 

0 

0 

0 

0 


453E-31 

0028366 

0031525 

0790777 

0796856 

0341976 

147E-39 


encounter 
angle (degrees ) 


Kl 


’ - tl 


T2 


0 

0.5221067 

66 .3312231 

12.8332741 

0.711E-35 

30 

2.4342451 

0.8039263 

7.4397717 

0.0017128 

60 

2 . 3829517 

0 . 8094406 

8 . 6272535 

0.0042339 

90 

2.0794163 

0.4690621 

16.9594727 

0.0442135 

120 

1.9938784 

0 . 2545378 

22.4892426 

0.  1164415 

150 

1.9387684 

0.2628353 

23.8461609 

0 . 0379958 

180 

0.7195223 

25.2119598 

14. 1219782 

0 . 312E-23 

encounter 

Encounter  Frequency  1.5  : 
Kl  Tl 

rads  per  sec 
T2 

J 

gle ( degrees ) 

0  0.5221067 

66 .3312231 

12.8332741 

0  .  143E-35 

30 

1.8784037 

0 .6553955 

5.4344263 

0 . 0000847 

60 

2.4894753 

0.5992758 

5  .  1160698 

0 . 0014724 

90 

2 . 5899792 

0 . 6663168 

3  .  1813316 

0 . 0241299 

120 

1.8784037 

0 .5395500 

10.4344263 

0 . 0028635 

150 

2.4478331 

0.5559916 

5.9208412 

0 . 0015984 

180 

0.7195223 

25.2119598 

14 . 1219782 

0 . 45  IE- 28 

table  vi 


Optimal  Controller  Parameters  for  Regular  Sea 

Sea  State  6 


Encounter  Frequency  0.2  rads  per  sec 


encounter 
angle (degrees ) 


.5221067 

0966187 

2665367 

3603411 

4170542 

4533644 

7195223 


66 . 3312231 
62.3766327 
55 .  1072540 
49 . 3086395 
46.4313202 
45.2830963 
25.2119598 


8332741 

3941193 

2263336 

7315979 

7203903 

6198730 

1219782 


Encounter  Frequencv  0.4  rads  per  sec 
encounter  K1  T1  T2 

angle (degrees ) 


5221067 

7536127 

1621914 

8681517 

6123600 

6361399 

7195223 


66  .  3312231 
12.4472111 
1.7206783 
28.5688019 
28.0108032 
27.8540497 
25 . 2119598 


Encounter  Frequencv  ( 
encounter  K.1  Tl 

angle (degrees ) 


rads 


5221067 

0449467 

0594482 

1537333 

1456242 

9752455 

.7195223 


66.3312231 

1.2405663 

0.1432046 

1.1794500 

0.9514294 

0.8255181 

25.2119598 


8332741 

8151121 

6173820 

3736725 

1387405 

7665482 

1219782 

per  sec 
T2 


8332741 

5909785 

5202637 

9987049 

3327799 

5351868 


,  322E-33 
622314 
84099 
270265 
004774 
396385 
678E-28 


.  122E-31 
3798 
143427 
7779746 
8634676 
2140626 
345E-25 


806E-35 

0113364 

0126108 

2766824 

3125796 

1364259 


14.1219782  0.112E-23 


encounter  K1 

angle (degrees ) 

0  0.5221067 

30  2.4413719 

60  2.4142313 

90  2.0592680 

120  2.0695038 

150  1.9513931 

180  0.7195223 

Encounter  Fi 
encounter  K1 

angle (degrees ) 

0  0.5221067 

30  1.7606564 

60  2.5002985 

90  2.6015043 

120  2.1809998 

150  2.4302263 

180  0.7195223 


Encounter  Frequency  0.75  rads 


per  sec 
T2 


66 . 3312231 
0 .7599964 
0 . 7818718 
0 . 3840635 
0 . 2843146 
1.2351496 
25.2119598 


12.8332741  0.753E-32 

7.3472633  0.0068484 

8.5500679  0.0169234 

19.2782745  0.1523162 

22.5692444  0.4640285 

23.5140228  0.1518951 

14.1219782  0.691E-28 


’requeue^ 


rads 


per  sec 
T2 


66  .  3312231 
0.4668925 
0.5747030 
0.6543975 
0.4698086 
0.5666089 
25.2119598 


12.8332741 
12.3475494 
5.3262844 
3 . 2408133 
10. 9160814 
6.0247307 
14. 1219782 


,  321E-33 
0003390 
0058868 
0946725 
0114672 
.0063947 
344E-23 
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TABLE  VII 


Optimal  Controller  Parameters  for  Pv.egular  Sea 

Sea  State  7 


Encounter  Frequency  0.2  rads  per  sec 


encounter 
angle (degrees ) 


0.5221067 
1.1839037 
1.3688898 
1.4492817 
1.4652090 
1.4653692 
0 .7195223 


3312231 

3722687 

5752258 

9667511 

2668457 

8378906 

2119598 


Encounter  Frequency 
encounter  K1  T ! 

angle (degrees ) 


rads 


0  0.5221067 

30  0.6452138 

60  2.2381306 

90  2.7780743 

120  2.5894566 

150  2.6306906 

180  0.7195223 

Encounter  Fi 
encounter  K1 

angle (degrees ) 

0  0.5221067 

30  2.0501146 

60  1.0624285 

90  2.1488247 

120  2.1970577 

150  2.0124264 

180  0.7195223 


3312231 

7243542 

0741718 

6558838 

1353302 

8760681 

2119598 


8332741 

3912964 

9693604 

6652069 

1322327 

1106033 

1219782 

per  sec 
T2 


8332741 

6571761 

7088852 

9282227 

0156784 

9550962 

1219782 


r  requenc^ 


rads 


66 . 3312231 
1.2125063 
0.  1251755 
0.9770966 
0.8763-103 
0.8064904 
25.2119598 


12.8332741 
5 . 5930214 
18.0873718 
8 . 1855469 
8.4206886 
8.9983368 
14.1219782 


encounter 
angle ( degrees ) 


Encounter 

K1 


Frequency 


9  HE-  35 

62723 

53816  ' 

31325 

58828 

036665 

912E-32 


341E-30 

.0524352 

22498 

95129 

098495 

6211755 

176E-35 


0.  157E-31 
0.0346353 
0.0386276 
0.6328694 
0.9178923 
0.4150548 
0 . 44  IE- 28 


0.75  rads  per  sec 


0 

0.5221067 

66.3312231 

12.8332741 

0 . 238E-28 

30 

2.4456072 

0.8701959 

7 . 6628313 

0 . 0209509 

60 

2.4188919 

0 . 8008499 

8.6116581 

0.0517269 

90 

1.9072247 

0 . 1301596 

28 .8665771 

0.4958954 

120 

2.2209492 

0 . 3453745 

22.8497772 

1.398130 

150 

2.0318069 

0.2577103 

23.6702576 

0.4641950 

180 

0 .7195223 

25 .2119598 

14.1219782 

0 . 542E-21 

Encounter  Frequency  1.5 

rads  per  sec 

ncounter 

K1 

T 1 

T2 

J 

le ( degrees  ) 

0 

0.5221067 

66  .  3312231 

12.8332741 

0.712E-32 

30 

2.0834208 

0.4371362 

15  .  1762695 

0 . 0010384 

60 

2.4635468 

0 . 5746776 

5 .3325920 

0.0180051 

90 

2 .6105270 

0 . 6474890 

3.4996614 

0 .2765892 

120 

2.  1780138 

0.4588630 

11.4343033 

0.0352355 

150 

2.4314880 

0 .5580992 

6.0874767 

0.0195939 

180 

0 .7195223 

25.2119598 

14. 1219782 

0.413E-26 

TABLE  VIII 


Optimal  Controller  Parameters  for  Regular  Sea 

Sea  State  9 


Encounter  Frequency  0.2  rads  per  sec 


encounter 
angle (degree 

K1 

s) 

Tl 

T2 

J 

0 

0.5221067 

66.3312231 

12.8332741 

0.493E-31 

30 

1.2344990 

97.2906342 

53.3980713 

31.88802 

60 

1.5019569 

77 . 3820190 

45 .6973572 

30 .87321 

90 

1.5020103 

77.2400513 

45 .6112671 

30.87321 

120 

1.5420017 

51.2350922 

23.8943634 

21. 75188 

150 

1.4879055 

44.4281464 

15.2255249 

8 . 599576 

180 

0.7195223 

25 .2119598 

14.1219782 

0.810E-39 

Encounter  F 

requency  0.4  ■ 

rads  per  sec 

encounter  K1  tl  T2  J 

angle (degrees ) 


0 

0.5221067 

66.3312231 

12.8332741 

0.743E-33 

30 

0.1618259 

95.8238471 

18 .3397064 

8 . 103665 

60 

2.3103380 

0 .9894915 

12.5459442 

20.30595 

90 

3.3482129 

2 .6129665 

4.7000313 

7.494904 

120 

2.5573973 

28 . 7854462 

12 . 1667938 

3 . 113852 

150 

2.6185656 

27 .9137726 

9 . 3553696 

1.324692 

180 

0.7195223 

25.2119598 

14.1219782 

0.534E-33 

Encounter  Fre 

quency  0 . 6 

rads  per  sec 

encounter  K1  tl  T2  J 


angle (degrees ) 


0 

0.5221067 

66.3312231 

12.8332741 

0.478E-31 

30 

2.0493793 

.  1.2093735 

5 . 6404839 

0.0820505 

60 

1.1003008 

0.1920759 

18 . 7369995 

0.0919840 

90 

2.0454016 

0 . 5407953 

15.4096680 

1.031631 

120 

2.2776527 

■  0.7766371 

10.4350891 

2.069717 

150 

2.0829115 

0.7725539 

9.8721237 

0.9771649 

180 

0.7195223 

25.2119598 

14.1219782 

0. 117E-29 

encounter 

Encounter  Frequency  0.75 
K1  Tl 

rads  per  sec 
T2 

J 

angle (degree 

!S) 

0 

0.5221067 

66 . 3312231 

12.8332741 

0.872E-35 

30 

2.4433956 

0 .8601446 

7  .  7009745 

0.0497640 

60 

2.4255047 

0 . 7809458 

8 . 7044067 

0. 1226490 

90 

2 . 2933350 

0.4751374 

16 .7149658 

0.7522082 

120 

2.3928595 

0.4103206 

22 . 9515076 

3 . 146294 

150 

2.1456175 

0.2891768 

23.5924225 

1.097565 

180 

0.7195223 

25 .2119598 

14.1219782 

0 . 6 11E-26 

encounter 

Encounter  F 
K1 

requenc^l .  5  : 

rads  per  sec 

T2 

J 

angle (degrees  ) 

0 

0.5221067 

66  .  3312231 

12.8332741 

0.242E-33 

30 

2.0808840 

0.4744592 

16  .  3423157 

0.0024732 

60 

2.4694090 

0.5775681 

5.4030972 

0 . 0042757 

90 

2 . 6412249 

0 .6233797 

3 . 9802380 

0.6127556 

120 

2 .  1955976 

0 . 4414446 

12 . 1491365 

0 . 0084519 

150 

2.4388838 

0 . 5490999 

6 . 1774960 

0 . 0046688 

180 

0.7195223 

25.2119598 

14.1219782 

0.420E-30 

30 


Figure  A.  2  Yaw  vs  Time,  Sea  State  A. 
counter  Frequency  1.5  rads  per  sec,  Encounter  Angle  120 


Rudder  vs  Time,  Sea  State  4. 

5  rads  per  sec,  Encounter  Angle  120 


ec 


Figure  4. 5  Rudder  vs  Time,  Sea  State  6. 
counter  Frequency  1 . b  rads  per  sec,  Encounter  Angle  120 


Encounter  Freque 


Figure  A.  9  Rudder  vs  Time,  Sea  State  9. 
Encounter  Frequency  l.b  rads  per  sec,  Encounter  Angle  120 


Cos 

rads 


re  A. 11  Cost  vs  Tl,Sea 
ency  1.5  rads  per  sec, 


Figure  4-12  Cost  vs  T2,Sea  State  4. 

Encounter  frequency  1.5  rads  per  sec,  Encounter  Angle  120 


Figure  4.13  Cost  vs  T2,Sea  State  6. 
requency  1.5  rads  per  sec.  Encounter  Angle  120 


Figure 
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V.  IRREGULAR  SEAS  -  CONTROLLER  DESIGN 


The  major  characteristic  of  the  sea  is  its  irregularity. 
This  irregularity  can  be  described  by  statistical  methods  by 
assuming  that  a  large  number  of  regular  (sinusoidal)  waves 
having  different  wavelengths,,  directions,  phases  and  ampli¬ 
tudes  are  superimposed  to  form  the  randomly  varying  sea. 

The  presence  of  the  irregular  sea  was  obtained  by 
coupling  a  sea  state  generator  program  to  the  FORTRAN 
program  as  is  indicated  in  Appendix  D.  The  sea  state  gener¬ 
ator  program  generates  added  mass  and  added  inertia  values 
as  function  of  the  encounter  frequency  and  also  calculates 
forces  and  moments  imparted  to  the  shiphull  by  the  sea.  The 
forces  and  moments  are  stored  in  a  look  up  table  which  was 
coupled  to  the  equations  of  motion.  The  irregular  sea  waves 
impinging  on  the  ship  contain  the  total  energy  density  spec¬ 
trum  composed  of  many  frequencies  and  the  ship  responds  to 
an  average  value  of  added  mass  and  added  inertia,  while  in 
the  regular  sea  the  added  mass  and  added  inertia  were  known 
for  a  given  encounter  frequency.  We  decided  to  use  values 
for  added  mass  and  added  inertia  corresponded  to  encounter 
frequency  0.75  rads  per  sec,  since  the  energy  density  is 
maximum  in  the  vicinity  of  this  frequency  [Ref.  5].  This 
frequency  gave  us  values  representative  of  an  average  value 
for  added  mass  and  added  inertia. 

The  controller  used  for  this  study  was  the  controller 
described  in  Chapter  4  (Figure  4.1).  The  optimized 
controller  parameters  and  the  cost  J  for  sea  states  4,  6,  7, 
9  and  0°,  30°,  60°,  90°,  120°,  150°,  180°  encounter  angles 
are  indicated  in  Table  IX. 

Studying  the  Table  IX  we  can  draw  the  following 
conclus ions : 
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•  For  sea  states  6,  7,  9,  the  higner  the  sea  state  the 

higher  the  cost,  for  every  particular  encounter  angle. 

•  Comparing  costs  for  sea  states  4  and  6  we  discover  some 
anomaly.  The  cost  for  a  specific  encounter  angle  in  sea 
state  4  is  higher  than  the  cost  for  the  same  encounter 
angle  in  sea  state  6.  Logically,  we  expect  higher  cost 
for  higher  sea  state. 

•  The  reason  for  this  anomaly  may  be  the  method  we  used 
in  order  to  obtain  the  added  mass  and  added  inertia 
values.  The  average,  we  consider,  might  not  represent 
the  actual  average. 

Appendix  E  provides  the  computer  program  necessary  to 
achieve  the  system's  response.  Some  typical  responses  are 
indicated  in  Figures  5.1,  5.2,  5.3,  5.4. 


TABLE  IX 

Optimal  Controller  Parameters  for  Random  Sea 


Sea  State  4 


encounter 

K1 

T 1 

T2 

J 

angle ( degrees 

) 

0 

0.6021814 

60 . 30849 

10 . 02579 

0 . 342E-24 

30 

1.5121580 

89 . 85324 

19 . 85960 

0 . 10719 

60 

0.6298036 

79.07199 

10.29221 

0 . 054196 

90 

0 . 6452737 

82 . 68692 

10 . 79342 

0 .076713 

120 

0 . 7485995 

85 . 77544 

12.37746 

0 . 137624 

150 

0 .9101038 

92 .71379 

15 .21078 

0 .012319 

180 

0.6021814 

60.30849 

10.02579 

0 . 189E-28 

Sea 

.  State  6 

encounter 

K1 

T 1 

T2 

J 

angle ( degrees 

) 

0 

0.6021814 

60 . 30849 

10.02579 

0 . 172E- 34 

30 

1.8743490 

61.82320 

32.22498 

0.044758 

60 

0 .8662014 

90.72922 

14.44058 

0 . 028238 

90 

0 . 7370305 

84.08502 

12 . 17457 

0.018541 

120 

2 . 6737600 

138 . 06650 

48 . 52447 

0 . 065028 

150 

0.4874309 

77 . 86977 

41.73848 

0 .012478 

180 

0.6021814 

60 . 30849 

10 . 02579 

0. 767E-23 

Sea 

State  7 

encounter 

K1 

T1 

T2 

j 

angle (degrees 

) 

0 

0. 6021814 

60 . 30849 

10 . 02579 

0 . 142E- 34 

30 

1.7232471 

66.44597 

27 . 32489 

0.41237 

60 

1.8508301 

91.39410 

36 .91092 

0 .377894 

90 

3.7642412 

62.27244 

86.58169 

0.258193 

120 

1.8482047 

99.64923 

91.28040 

0.225806 

150 

0.8519831 

67.02328 

64.96774 

0 . 037724 

180 

0.6021814 

60.30849 

10.02579 

0.811E-21 

Sea 

State  9 

encounter 

x  K1 

T1 

T2 

j 

angle ( degrees 

) 

0 

0.6021814 

60.30849 

10 . 02579 

0 . 711E-35 

30 

3 . 1908160 

138 . 56101 

71.44171 

1.  306741 

6  0 

3.0888780 

152  -  01630 

72 .66160 

0 . 961709 

90 

3.2440700 

121. 13566 

105.98480 

0.523769 

120 

1.5461040 

111.49130 

99 . 64659 

0 .2101076 

150 

0.3758357 

73 . 88629 

35  .  17305 

0 . 5007348 

180 

0.6021814 

60 . 30849 

10.02579 

0. 12  IE  -  3  1 
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Figure  5.2  Rudder  vs  Time,  Sea  State 
Encounter  Angle  60°. 


Figure  5.3  Yaw  vs  Time, 
Encounter  Angle 
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VI.  MINIMIZATION  SUBROUTINE  FOR  ONBOaRO  USE 

A.  GENERAL 

As  is  mentioned  earlier  in  Chapter  1  the  function  mini¬ 
mization  subroutine  used  for  these  studies  was  BOXPLX.  This 
subroutine  will  find  the  minimum  of  any  function,  linear  or 
nonlinear,  subject  to  explicit  constraints  of  the  variables 
or  implicit  constraints  on  functions  of  the  variables.  It 
will  handle  a  maximum  of  25  variables  but  can  handle  up  to 
50  variables  with  user  modification. 

The  variables  in  BOXPLX  are  allowed  to  move  within  a 
feasible  region  ( n- dimens ional  space,  where  n  is  the  number 
of  variables)  defined  by  upper  and  lower  bounds  on  their 
values.  The  choices  •  for  upper  and  lower  bounds  for  the 
parameters  are  based  on  an  understanding  of  the  function  of 
each  coefficient  of  the  system.  Experience  indicates  that 
while  accurate  selections  of  these  bounds  are  not  necessary, 
intelligent  selection  of  these  as  well  as  the  starting 
values  (guesses)  can  considerably  reduce  the  computer  number 
of  trials  needed  for  solution  convergence.  This  conclusion 
was  drawn  trying  to  obtain  the  controller  parameters  for 
Tables  V,  VI,  VII,  VIII  in  Chapter  4.  The  function  minimiza¬ 
tion  subroutine,  when  starting  the  minimization  process  with 
arbitrary  chosen  guesses,  required  more  than  100  trials  for 
convergence  while  by  choosing  guesses  close  to  the  optimal 
parameter^  required  more  than  50  and  less  than  100  trials. 
Considering  that  every  trial  lasted  600  seconds  (10  minutes) 
and  the  function  minimization  subroutine  requires  about  60 
samples  (trials)  before  telling  us  it  had  found  the  minimum 
this  would  mean  10  hours  for  the  control  to  adjust  itself. 
For  obvious  reasons,  such  operation  is  not  acceptable  for  on 
board  use . 


For  obvious  reasons,  such  operation  is  not  acceptable  for  on 
board  use. 

B.  ATTACKING  THE  PROBLEM 

We  started  to  investigate  ways  to  improve  this.  These 
efforts  include: 

•  Finding  a  more  efficient  function  minimization 
subroutine 

•  Studying  the  flatness  of  the  cost  surface 

•  Reducing  sampling  time 

C.  SOLVING  THE  PROBLEM 

Switching  to  another  function  minimization  subroutine  we 
found  that  the  new  one  (ZXMWD)  suffered  from  the  same  disad¬ 
vantages  . 

The  experiments  carried  out,  more  than  two  hundred, 
indicated  that  the  cost  surface  is  really  flat.  The  BOXPLX 
after  a  few  trials  started  to  focus  on  the  minimum  but 
before  it  converged,  it  needed  more  than  50  trials,  even  if 
the  guesses  were  close  to  the  optimal.  The  reason  is  the  way 
BOXPLX  itself  tries  to  find  the  minimum  of  a  function  of  NV 
variables.  It  converges  when  the  cost  FE  remains  unchanged 
for  2"NV  consecutive  trials  with  accuracy  10s.  An  effort  to 
modify  this  termination  criterion  in  terms  of  the  consecu¬ 
tive  trials  was  successful.  Table  X  indicates  the  comparison 
between  modified  and  unmodified  BOXPLX,  for  sea  state  6, 
encounter  frequency  1.5  rads  per  sec  and  encounter  angle 
120°.  As  we  can  observe  in  Table  X  the  cost  in  each  case 
remains  almost  the  same  while  the  trials  required  for 
convergence  are  dependent  on  the  guesses  made  and  the  termi¬ 
nation  criterion  established. 

The  value  of  the  cost  is  in  general  the  summation  of 
incremental  contributions  for  each  integration  step  and  is 


TABLE  X 

First  Modification  in  30XPLX 


BOXPLX 

Guesses 

Trials 

Terminat ion 
Criterion 

Cost 

Unmodified 

Arbitrary 

100 

6 

0.01146720 

Unmodified 

Optimal 

72 

6 

0 .01146720 

Modified 

Arbitrary 

42 

3 

0 .01146812 

Modified 

Optimal 

11 

3 

0.01146720 

Modified 

Arbitrary 

7 

2 

0.01157926 

Modified 

Optimal 

1 

2 

0 .01146720 

therefore 

dependent  on 

the 

total  time  of 

the  simulation 

This  is  important  in  that  the  optimal  gain  coefficients 
arrived  at  in  this  manner  are  not  optimal  for  steady-state 
performance  but  only  for  the  time  frame  covered.  This  should 
be  adequate  provided  the  time  frame  selected  is  long  rela¬ 
tive  to  the  time  required  for  the  initial  condition  response 
to  die  out.  This  is  the  reason  Reid  has  chosen  time  frame 
600  seconds  [Ref.  1,2].  Of  course  this  time  period  is  large 
and  we  expect  steady-state  behaviour  faster  than  10  minutes. 
Simulation  studies  indicate  that  the  ship,  controlled  by  the 
controller  described  in  Figure  4.1,  reaches  the  steady-state 
situation  in  less  than  100  seconds.  So,  we  can  reduce  the 
600  seconds  time  frame  to  200  seconds,  safely.  This  is  very 
important  since  now  the  modified  function  minimization 
subroutine  BOXPLX,  uses  samples  of  200  seconds  long  instead 
of  600  seconds,  converges  in  less  than  30  minutes  which  is 
reasonable  for  on  board  use.  Since  the  value  of  the  cost  is 
dependent  on  the  time  frame  taken  we  expect  reduced  cost  for 
200  seconds  samples  long,  in  any  case. 

As  we  discussed  earlier  the  BOXPLX  compares  the  consecu¬ 
tive  trials  with  accuracy  10s.  Since  we  do  not  need  so  big 
accuracy  a  second  modification  is  necessary.  We  decided  to 
change  the  existing  accuracy  to  10“. 

Table  XI  indicates  the  trials  required  for  BOXPLX 
convergence,  after  the  second  and  last  modification,  for  sea 
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TABLE  XI 

Second  Modification  in  BOXPLX 


Guesses 


Trials  Termination 

criterion 


Arbitrary  17  3 
Optimal  8  3 
Arbitrary  5  2 
Optimal  1  2 


Cost 


0.003999837 

0.003995277 

0.004024354 

0.003995247 


state  6,  encounter  frequency  1.5  rads  per  sec  and  encounter 
angle  120°.  By  comparing  Tables  X  and  XI  we  can  see  the 
further  improvement  for  the  function  minimization  subroutine 
convergence.  The  difference  in  the  cost  is  due  to  the 
different  time  frames  used.  The  modified  function  minimiza¬ 
tion  subroutine  BOXPLX  is  indicated  in  Appendix  F. 


VII.  ADAPTIVE  CONTROL 


A.  NECESSITY  OF  ADAPTIVITY 

The  plant,  the  system  which  is  supposed  to  be 
controlled,  is  normally  exposed  to  a  time  varying  environ¬ 
ment.  The  ship's  speed,  the  wave  encounter  angle  ~nd  the  sea 
state  are  changing  drastically  during  the  seaway.  Since 
changes  encountered  are  not  completely  predictable  an 
optimum  preprogrammed  time-varying  controller  is  not 
possible . 

If  we  assume-and  it  is  apparent  from  the  previous 
discussion  in  this  study-that  no  feasible  fixed  parameter 
controller  provides  acceptable  response  over  the  entire 
performance  spectrum,  it  is  necessary  that  some  means  is 
provided  for  adjusting  controller  parameters  according  to 
the  sea  conditions  and  ship’s  operational  characteristics. 
Adaptive  control  is  thus  an  effort  to  extend  basic  optimum 
control  concepts  to  these  studies. 

B.  CANDIDATE  ADAPTIVE  SCHEMES 

Since  we  are  looking  for  1%  or  2%  savings  in  fuel  cost , 
we  have  to  feed  the  system  with  precise  information.  Exact 
knowledge  of  the  sea  state,  the  wave  encounter  angle  and  the 
ship's  ground  speed  is  vital  for  this  purpose.  Currently  the 
Navy  is  involved  in  a  program  that  will  provide  precision 
navigation  data.  Garcia  [Ref.  5],  provides  very  good  infor¬ 
mation  on  this  subject. 

An  adaptive  control  scheme  is  indicated  in  Figure  7.1. 
Once  the  adaptive  part  of  the  scheme  is  set  up  the  sea 
state,  encounter  wave  angle  and  ship’s  speed  are  fed  to  the 
filters  box  by  the  appropriate  sensors.  The  filters  box 


includes  predetermined  optimal  filter  sets  for  different 
discrete  sea  states,  wave  angles  and  ship's  speed.  Actually, 
it  is  a  look  up  table  similar  to  chose  of  Tables  V,  VI, VII 
and  VIII.  The  output  of  the  filter's  box  is  a  filter  which 
corresponds  to  discrete  conditions  close  to  those  fed  by  the 
sensors.  The  function  minimization  subroutine  accepting  the 
rudder  angle,  yaw  error  and  the  predetermined  filter  set  as 
initial  guesses  tries  to  obtain  the  optimal  filter  for  the 
exact  sea  state,  wave  encounter  angle  and  ship's  speed.  At 
the  same  time  the  plant  is  controlled  by  the  controller  with 
the  optimal  predetermined  set  of  parameters  for  conditions 
close  to  the  actual.  When  the  function  minimization  subrou¬ 
tine  reaches  the  minimum,  it  supplies  the  controller  with  a 
new  set  of  parameters  which  is  the  optimal  set  for  the 
present  conditions. 

But,  what  happens  if  either  some  or  all  the  sensors 
provide  new  inputs  to  the  system?  A  decision  device  placed 
after  the  sensors  decides  whether  or  not  the  change  is 
appreciable.  This  device  compares  the  current  conditions 
with  those  used  to  obtain  the  controller  which  currently 
governs  the  system.  If  the  change  is  higher  than  some 
desired  percentage  then  a  new  predetermined  filter  is  passed 
to  the  controller  and  the  function  minimization  subroutine 
tries  to  find  the  optimal  controller  parameters  for  the  new 
situation. 

From  the  scheme  of  Figure  7.1  we  can  eliminate  the 
predetermined  filters  device  which  provides  a  less  expensive 
system.  In  this  case  the  function  minimization  subroutine 
will  need  a  much  longer  time  to  determine  the  optimal  filter 
for  rapid  changes  in  course  and  speed,  even  if  we  assume 
that  rapid  changes  in  sea  state  don't  occur.  Finally,  the 
function  minimization  subroutine  might  work  continually  "on 
line"  as  is  indicated  in  Figure  7.2.  In  this  scheme  a  new 
controller  set  is  obtained  in  every  subroutine  iteration  and 


some  constraints  for  the  controller  parameters  are  necessary 
in  order  to  avoid  operation  in  unstable  regions.  Thus,  the 
filter  supplied  by  the  subroutine  for  controlling  the  plant 
must  be  tested  for  characteristic  equation  roots  in  the 
right  half  S-plane.  In  that  case  a  default  option  must 
exist,  and  a  special  device  for  such  purposes  is  necessary. 
This  device  will  permit  change  in  the  filter  parameters  only 
when  the  new  set  still  keeps  the  system  in  a  stable  situ* 
at  ion . 

Whichever  adaptive  scheme  we  adopt,  we  must  provide  for 
manual  operation  for  the  system.  Manual  operation  will  be 
desired  in  the  following  cases: 

•  Arriving  ports 

•  Leaving  ports 

•  Restricted  waters 

•  Avoid  collision  in  open  seas 

•  Computer  down 

For  the  first  two  situations,  since  we  usually  expect  no 
heavy  seas,  the  optimal  filter  for  sea  state  1  seems  to  be 
more  suitable.  Also,  this  filter  can  serve  as  initial  condi* 
tion  for  the  adaptive  schemes  described  before,  when  we  are 
leaving  ports.  For  the  rest  of  the  situations  the  last 
operating  optimal  filter  is  the  most  appropriate. 


VIII.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

The  principal  conclusions  from  this  study  of  the  SL-7 
containership  as  they  related  to  steering  may  be  stated  as 
follows : 

•  It  is  evident  that  a  control  system  which  provides  the 
ship  heading  and  simultaneously  reduces  the  propulsive 
losses  does  exist  and  therefore  such  a  controller  saves 
fuel.  We  can't  conclude  how  much  the  savings  are, 
since  there  is  no  reference  for  comparison  between  the 
conventional  autopilot  and  che  autopilot  which,  in 
addition,  holds  the  potential  for  reducing  the  propul¬ 
sive  losses.  The  literature  says  that  savings  1%  or  2Z 
is  possible. 

•  An  adaptive  controller,  that  minimizes  propulsion 
losses  as  ship  characteristics  and  environmental  condi¬ 
tions  change,  may  be  designed  using  a  self-optimalizing 
technique  employing  a  suitable  performance  criterion. 

•  Studying  every  particular  situation  we  conclude  that 
the  cost  surface  is  flat  and  therefore  accurate  deter¬ 
mination  of  the  controller  parameters  is  not  required. 
So,  if  we  decide  to  use  the  adaptive  control  scheme  of 
Figure  7.1  we  don't  have  to  store  every  particular 
filter  in  the  look  up  table  since  one  filter  may  be 
suitable  for  different  ship  characteristics  and  sea 
conditions . 

•  The  weighting  factor  X  used  in  the  performance 
criterion  equation  3.5  plays  an  exceptional  and  impor¬ 
tant  role  in  the  optimal  controller  parameters  determi¬ 
nation.  However,  it  is  obtained  from  studies  based  on 
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many  assumptions  and  therefore  isn't  predicted 
rately.  As  a  consequence  the  controller  found  does  not 
minimise  added  res  is  t ance  unless  the  proper  value  of 
the  weighting  factor  X  has  been  used. 

•  The  method  used  to  obtain  the  average  for  the  added 
mass  and  added  inertia  for  the  irregular  seas  studies 
might  not  represent  the  actual  average. 

•  The  function  minimization  subroutine  BOXPLX  after  two 
modifications  seems  to  be  pretty  suitable  for  on  board 
use  working  as  a  main  part  of  the  adaptive  scheme. 

B.  RECOMMENDATIONS  FOR  FUTURE  STUDIES 

The  following  recommendations  for  future  work  to  gain  a 
fuller  and  deeper  understanding  of  the  problem  are  made  as 
follows : 

•  Some  studies  are  necessary  in  order  to  investigate  why 
part  of  the  obtained  controllers  in  Tables  V,  VI,  VII, 
VIII  are  lag  and  part  of  them  are  lead  filters. 

•  As  we  stated  in  chapter  4  the  yaw  and  rudder  excursions 
in  Figures  4.2  through  4.9  are  less  than  1*.  It  is 
necessary  to  investigate  the  reason  for  that.  It  might 
be  because  of  the  optimization  of  the  filter  or  the 
forces  and  moments  have  not  been  sealed  properly. 

•  It  is  necessary  to  find  out  the  appropriate  average 
values  for  the  added  mass  and  added  inertia  before  we 
attempt  further  studies  in  irregular  seas. 

•  The  full  hydrodynamic  coefficients  for  the  SL-7  are 
necessary  in  order  to  develop  and  include  the  surge 
equation  in  our  ship  model.  So  far  we  ignore  the  surge 
equation  and  we  assumed  constant  ship  speed  while  in 
reality  the  added  resistance  due  to  steering  must 
reduce  the  ship  speed. 
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•  By  developing  the  surge  equation  in  our  ship  model  we 
may  be  able  to  determine  a  good  value  for  the  weighting 
factor  to  be  used  in  the  performance  criterion. 

•  Also,  with  the  surge  equation  available  in  our  model  we 
should  be  able  to  calculate  actual  energy  losses  and 
savings  in  fuel. 

•  As  we  stated  in  chapter  6  the  time  frames  (samples) 
used  by  the  BOXPLX  must  be  long  relative  to  the  time 
required  for  the  initial  condition  response  to  die  out. 
Reid  [Ref.  1,2]  has  chosen  time  frame  600  seconds  long. 
This  time  frame  is  long  and  it  is  necessary  to  find  out 
the  sufficient  time  frame  since  the  time  required  by 

.  the  function  minimization  subroutine  for  convergence  is 
directly  proportional  to  the  sample  duration. 
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APPENDIX  A 

NOMOTO  THIRD  ORDER  MODEL  DETERMINATION 

//NOMOTO  JOB  (XXXX.XXXX) , 'RESEARCH' ,CLASS=J 

//-•MAIN  ORG  =  NPGVMI .  XXXXP 

//  EXEC  FORTXCG , PARM . FORT= ' OPT ( 2 ) ’ , IMSL=DP , REGION = 10 24K 

//FORT. SYS IN  DD  * 

C  IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE  BEEN 

C  OBTAINED  CHANGE  XS(*)  TO  X(*)  AND  DELETE  XU  ( •• )  .  AND  XL(*) 

DIMENSION  XS (4 ) , XU (4 ) , XL (4 ) 

XS(I)=0. 1 
XS(2)=I5 . 13 
XS ( 3 ) = 15 . 675 
XS(4)=9 .014 

C  XS(I)  IS  THE  STARTING  GUESS 

C  XL ( I )  IS  THE  LOWER  LIMIT  FOR  THE  I ’ TH  VARIABLE 

C  XU (I)  IS  THE  UPPER  LIMIT  FOR  THE  I ' TH  VARIABLE 

XL(1)= .01 
XU ( 1 ) = 1 . 0 
XL ( 2 ) = 1 . 0 
XU ( 2 )  =  20 . 

XL ( 3 ) = 1 . 0 
XU ( 3 ) = 100 . 

XL (4) =1.0 
XU (4) =100.0 

C  A  DESCRIPTION  OF  THE  FOLLOWING  PARAMETERS 

C  IS  DISCUSSED  IN  BOXPLX 
R=9  . / 13  . 

NTA= 1000 
NPR=  0 
NAV  =  0 


THE  FOLLOWING  STATEMENT  MUST  BE  CHANGED  TO 
CALL  PLANT (XX) 

IF  ONLY  SIMULATION  IS  WANTED 

CALL  BOXPLX ( NV , NAV , NPR , NT A , R , X S , I P , XU , XL , YMN , I ER ) 
WRITE  (6,25) 

j  FORMAT (IX,*  OPTIMAL  GAINS',/) 

DO  30  1=1,4 
)  WRITE (6,40)I,XS(I) 

)  FORMAT (IX,  ’ X( ’  ,12,  '  )='  ,  F14 . 7 ) 

WRITE  (6,77)  TDIFF 
7  FORMAT (IX, ’COSTS' ,F14. 7) 

STOP 


SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 
SUBROUTINE  PLANT (XX) 

COMMON  TDIFF 

REAL “8  L , L2 , L3 , L4 , L5 , L6 , RXR , RXI , RYR , RYI , MZR , MZI , RX , RY 
REAL " 8  X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT , T X , TY 
REAL  ""8  TIME , ETIME , XUDOT , XU , XUU , XVR , XVV , XDD , WA , WE , RZ 
REAL" 8  YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT , TZ 
REAL" 8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL "8  RHO , IZ , FX , F Y , MZ , XP , MAS S , DELT 
REAL "8  YAWE , YAWC , D2 , D 

REAL" 8  K , TP  1 , TP2 , Z , XI ,X2 , X3 , X4 , X5 , YAW2 , S 
DIMENSION  XX (4) 

INITIAL  CONDITIONS  FOR  INTEGRATION 
SIMULATION  END  TIME  IN  SECONDS 
ETIME=  600 . 

TIME=0 .0 
ICOUNT  = 1 

INITIALIZE  THE  COST  FUNCTION 
TDIFF  =  0.0 

GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 
K=XX( I) 

Z=XX ( 2 ) 


TP1=XX(3) 

TP2=XX (4 ) 

C  X , XDOT , Y , YDOT  ARE  FIX  COORDINATES  ON  EARTH 
X  =  0.0 
Y=0. 0 
XDOT  =  0 . 0 
YDOT  =  0 . 0 

C  U , UDOT , V , VDOT  ARE  FIX  COORDINATES  ON  SHIP 
V  =  0.0 
UDOT  =  0 . 0 
VDOT  =  0 . 0 
YAW  =  0 . 0 
R=0 .0 
RD0T  =  0 . 0 

C  ORDERED  SPEED  IN  FEET/SEC 

C  38.81  FT / SEC=  23  KNOTS 
UC=38.81 

C  AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 
U  =  UC 

C  D  =  RUDDER  ANGLE 
D  =  0 . 0 
L-880.5 
L2=L**2 
L3  =L*L*L 
L4=L*L3 
L5=L*L4 
L6  =  L*L5 

C  SEA  DISTURBANCE 

C  FORCES  IN  X , Y  DIRECTION  COMPUTED  IN  FORCES 

C  MOMENTS  IN  Z 
FX  =  0 . 0 
FY  =  0 . 0 
MZ  =  0 . 0 

C  RXR=- . 15744D-05 

C  RXI-- . 1995 OD *06 
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c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


RYR=0 . 52365D-04 
RYI = 0 . 18699D+06 
MZR= - . 29870D+08 
MZI=- . 37751D+07 
RXR= - . 50230D-04 
RXI=0. 12712D+05 
RYR=  0 . 35290D+04 
RYI= - . 3L909D+05 
MZR=0.38826D+07 
MZI= - . 643 13D+  07 
RXR=  0 . 28540D+04 
RXI= - . 99574D+04 
RYR= - . 85441D+04 
RYI=0 .39595D+05 
MZR=- . 13014D+08 
MZI=0. 11348D+08 
RXR= - . 7  5  642D+  04 
RXI=0 .8349 7D +04 
RYR=0 .23379D+05 
RYI = - . 8 1502D*05 
MZR=0 .28622D+07 
MZI= - . 19388D+08 
RXR= - .379L6D+04 
RXI  =  0 . 16  38 1D  +  04 
RYR= - . 7  6647D  +  05 
RYI= - . 37685D+05 
MZR= - . 83915D+07 
MZ I = - .53176D+07 
RX=DSQRT(RXR**2+RXI**2 ) 

RY  =  DSQRT  ( RYR"  ••  2  *  RYI  *  *2  ) 

RZ  =  DSQRT  ( MZR'-V"2  +  MZ I  **2  ) 
TX=DATAN2(RXI,RXR) 
TY=DATAN2 (RYI , RYR ) 
TZ=DATAN2 (MZI ,MZR) 
SIGNIFICANT  WAVE  HEIGHT;  SEA 


STATE  1-0.32,2-0.75,3-2.5 
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C  4-5.0,5-7.0,6-10.0,7-17.5,8-20.5,9-27.0 
WA= 10 . 0 

C  ENCOUNTER  FREQUENCY  .1,. 2, .3, .4, .5, .6, .75, 1.0, 1.5, 2. 5 
WE  =  0 .4 

C  HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  AS  PARAMETERS 
RHO  =1.9876 

MASS  =  (  .  0044)"  (  •  5 ••RHO" L 3  ) 

I Z  = ( 0 . 0  00  2  8) * ( . 5 " RHO " L5 ) 

YAWE  =  0 . 0 
X1=0.0 
X2  =  0 . 0 
X3  =  0 . 0 
X4  =  0 . 0 
X5  =  0 . 0 
YAW2=0 . 0 
200  CONTINUE 

S=DSQRT (U**2  +  V" " 2 ) 

C  INPUT  YAW  COMMAND 
YAWC  =  0 . 0 

IF  (TIME. GE. 0.0)  YAWC= (1 . 0/57 . 296  ) 

C  ERROR  SIGNAL  TO  DRIVE  RUDDER  (YAW  ACTUAL  -  YAW  COMMAND) 
C  FOR  EQUATIONS  OF  MOTION. 

YAWE= YAW  -  YAWC 
D=YAWE 
C 

C  NOMOTO  3RD  ORDER  PLANT 
C 

C  ERROR  SIGNAL  TO  DRIVE  RUDDER  (YAW  COMMAND  -  YAW  ACTUAL) 
C  FOR  NOMOTO  MODEL. 

D2=YAWC- YAW2 
Xl= ( D2 - X2 ) / TP  1 
X3  =  K " ( Z*X 1 *  X2 ) 

X4= (X3-X5 ) /TP2 

C  AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 

XUDOT  =  (  -  .  0001 )  ■•'•'  (  .  5  "RHO" L 3  ) 
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XUU  =( -0.0003  )'-•(. .5  *RHO*L2  ) 

XU= ( -  0 . 025  3  ) * ( 0 . 5*RHO*L2*S ) 

XVR=  ( 0 . 00 3  9  )'•■(.  5  "RHO'vL3  ) 

XVV=  (  -  .  00 12 ) *  (  .  5  "RHO',fL2  ) 

XDD= ( -0 . 0005 )*( . 5  * R H 0 “ L 2  * S  * * 2 ) 

C  LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 

C  YV=(-0. 00758 )*( ,5*RHO*L2*S) 

YR=  (0 . 0023  )*(’.  5*RHO*L3*S  ) 

YD  =  ( 0 . 0  0 1 4  5  )  *  (  .  5  *  RHO  *L  2  *  S  **  2  ) 

YVVR= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

YVRR= (-0.008) *(.5 *RH0*L4 / S ) 

YWV  =  (  -  0 . 03  )*(  .  5  "RHO"L2 /  S  ) 

YRRR=  (0 . 003  )*(  .5*RHO*L5/S) 

YDDD= ( -0 . 0005 )*( . 5*RHO*L2*S**2 ) 

YVDOT  =-0.30908D+07 
YV =-0.8127 ID +04 
C  YVDOT  =  - . 36185D  +  07 

C  YV=- . 24757D+06 

C  YVDOT=- . 32890D+07 

C  YV  =- . 1177  5D  +  07 

C  YVDOT = -. 23038 D+ 07 

C  YV= - . 18267D+07 

C  YVDOT=- .59800D+06 

C  YV=- . 13260D+07 

C  MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV=  (  -0 . 00213  )'•'(  .  5*RHO*L3*S  ) 

G  NR= ( -  0 . 00 10  5 ) * ( . 5*RHO*L4*S ) 

ND= ( -0 . 0007 )*( . 5 “RHO *L  3  * S  ** 2 ) 

NVVR= ( -0 . 015 )*( . 5*RH0*L4/S ) 

NVRR=  (  -  0 . 00 8  )'••(.  5  -RHO -L5  /  S  ) 

NVVV= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

NRRR= (-0.006 )*( .5*RHO*L6/S) 

NDDD=  (0 . 0001) *(  .  5 *RHO *L3 * S -- * 2 ) 

C  NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE  CHANGED 
C  FOR  DIFFERENT  ENCOUNTER  ANGLE , SPEED , ENCOUNTER  FREQUENCY 
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C  NRDOT  = ( -  0 . 00  0  2  7  ) * ( . 5 *RHO*L5 ) 

C  SPEED=23  KNOTS , ENCOUNTER  ANGLE  =  60, ENCOUNTER  FREQ=0.75 
C  NRDOT=- .26251D-12 

C  NR=- . 53637D-09 

NRDOT = - .20125D-12 
NR=- . 94970D* 10 
C  NRDOT=- . 18671D- 12 

C  NR=- .46860D-11 

C  NRDOT = - . 14518D+12 

C  NR= - . 8  75  38D+ 11 

C  NRDOT = -.37261D+11 

C  NR= - . 69856D+ 11 

C  SETS  SEA  STATE  TO  ZERO 
FX  =  0  . 

FY  =  0  . 

MZ  =  0  . 

FX= WA*RX*DCOS ( WE*TIME+TX ) 

F Y= WA*RY*DCO S ( WE “TIME + TY ) 

MZ  =  W A " RZ " D  C  0  S (WE*TIME*TZ ) 

C  U  ACTUAL  SPEED 
C  UC  COMMANDED  SPEED 
C  XP  =  PROPELLER  THRUST 
X  P  =  -  XUU  -  U  C  *  *  2 
C  EQUATIONS  OF  MOTION 

C  UDOT=  (  ( XVR  +  MASS)“V"R  *  XUU-U  —  2  +  XW*V**2 
C  1  +  XDD“D“D  *  FX  +  XP  ) / (MASS-XUDOT ) 

VDOT=  ( YV*V  +  ( YR  -  MA  S  S  *U )  *R  +  YD“D  +  YWR*V**2*R 

1  +  Y  V  R  R  “  V  “  R  *  *  2  *  YVVV“'V"V"  “  3 

2  ♦  YRRR“R““3  +  YDDD“D““'3  +  FY  )  /  (MASS -  YVDOT ) 

RDOT= (NV*V  +  NR“R  +  ND*D  +  NVVR“V"“2“R  +  NVRR*V*R**2 
1  +NVW*V**3  *  NRRR--R““3  +  NDDD*D**3  ♦  MZ  )  /  ( IZ  -  NRDOT) 
C  WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.il)  GO  TO  50 
GO  TO  300 


C  CONVERT  RADIANS  TO  DEGREES 
50  YAWDEG=  YAW-57.296 
RDEG=R- 57 . 296 
RDDEG=RDOT*57 . 296 
DDEG=D'-57 . 296 
YAWC=YAWC*57 .296 
ICOUNT  =  1 

C  TEST  IF  WANT  TO  STOP 
300  IF  (TIME . GE . ETIME )  GO  TO  400 
C  INTEGRATION  STEP  SIZE  DELT 
DELT= 1 . 0 
C  INTEGRATION 

X2  =  X2  +X1 -DELT 
X  5  =  X5  +  X4 -DELT 
YAW2= YAW2+X5 -DELT 
U=U+UDOT*DELT 
V  =  V+ VDOT -DELT 
R=R+RDOT*DELT 
YAW = YAW + R'- DELT 

C  CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 
XDOT=U-DCOS  (YAW )  - V'-D-SIN (YAW ) 

YDOT=U*DS IN ( YAW ) ♦ V - DCOS ( YAW ) 

X=X*XDOT-DELT 
Y= Y  +  YDOT '-DELT 
TIME = TIME + DELT 
I COUNT = ICOUNT  +  I 
C  COST  FUNCTION 

TDIFF=TDIFF (  YAW-YAW2  )**2 
GO  TO  200 
400  CONTINUE 

C  WRITE (6,500)  TDIFF , K , Z , TP  I , TP2 
500  FORMAT (’  ’ , IX , ’  COST  = ’ , FI2 . 7 , 2X , ’  K  =’,FI0.7, 
I  ’  Z  = ’ , FI5 . 7 , ’  TP  1  =’  ,FI5.7,2X, ’  TP2  =’,F15.7) 
RETURN 
END 


72 


APPENDIX  B 


REGULAR  SEASTATE  FORMULATION 

/ / REGUGAINS  JOB  (XXXX , XXXX ),' RESEARCH CLASS= C 
//-MAIN  ORG=NPGVMl . XXXXP 

//  EXEC  FORTXCG , PARM. FORT= ' OPT (2 ) ' , IMSL=DP , REGION = 102 4K 
/ / FORT . SYSIN  DD  * 

C  IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE  BEEN 

C  OBTAINED  CHANGE  XS(*)  TO  X(*)  AND  DELETE  XU (*),  AND  XL ( •-  ) 

DIMENSION  XS ( 3 ) , XU ( 3 ) , XL ( 3 ) 

XS(1)=0 . 9650610 
XS(2)=0. 4500911 
XS(3)=5 .6194260 

C  XS(I)  IS  THE  STARTING  GUESS 

C  XL(I)  IS  THE  LOWER  LIMIT  FOR  THE  I ' TH  VARIABLE 

C  XU(I)  IS  THE  UPPER  LIMIT  FOR  THE  I’TH  VARIABLE 

XL ( 1 )  =  0 . 1 
XU  ( 1 )  =  4 . 0 
XL ( 2 )  =  0 . 1 
XU(2 )  =  15 . 0 
XL ( 3 ) = 1 . 0 
XU ( 3 )  =  2  5 . 0 

C  A  DESCRIPTION  OF  THE  FOLLOWING  PARAMETERS 
C  IS  DISCUSSED  IN  BOXPLX 
R=  9 . / 13  . 

NTA= 1000 
NPR= 100 
NAV  =  0 
NV=  3 
IP  =  0 

C  THE  FOLLOWING  STATEMENT  MUST  BE  CHANGED  TO 
C  CALL  PLANT ( X ) 

C  IF  ONLY  SIMULATION  IS  WANTED 


CALL  BOXPLX ( NV , NAV , NPR , NTA , R , XS , IP , XU , XL , YMN , IER ) 
WRITE  (6,25) 

25  FORMAT ( IX , ’  OPTIMAL  GAINS '  ,  /  ) 

DO  30  1=1,3 

30  WRITE (6,40)1, XS ( I ) 

40  FORMAT (IX, *  X ( ’  ,12,')=' , F14 . 7  ) 

STOP 

END 

SUBROUTINE  PLANT (XX) 

C  SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 
COMMON  TDIFF 
REAL'" 8  L,L2  ,L3  ,L4  ,L5  ,L6 

REAL "8  X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 
REAL"  8  TIME  ,  ETIME  ,  XUDOT  ,  XUU  ,  XVR ,  XW , XDD 
REAL -8  YV , YR , YD ,  YWR  ,  YVRR  ,  YVW  ,  YRRR ,  YDDD  ,  YVDOT 
REAL "8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL'"  8  RHO  ,  I Z  ,  FX ,  FY  ,  MZ  ,  XP  ,  MAS  S  ,  DELT  ,  MZ I ,  WA  ,  WE 
REAL" 8  DYAWE , YAWE , YAWC , ISE , ISR , LAMDA , D , RYR , RYI , MZR 
REAL'"  8  K1 ,  T1 ,  T2  ,  D  ,  X2  ,  DX2  ,  S  ,  RX  ,  RY  ,  RZ  ,  TX  ,  TY  ,  TZ  ,  RXR  ,  RXI 
DIMENSION  XX ( 3 ) 

C 

C  CLOSE  LOOP  ANALYSIS  WITH  FILTER 

C 

C  INITIAL  CONDITIONS  FOR  INTEGRATION 

C  SIMULATION  END  TIME  IN  SECONDS 
ETIME  =  600 . 0 
TIME=0 .0 
ICOUNT= 1 

C  INITIALIZE  THE  COST  FUNCTION 
ISE=0 . 0 
ISR=0 . 0 
TDIFF  =  0.0 
LAMDA=  8 . 128 

C  GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 
K1=XX ( 1 ) 
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T1  =  XX ( 2 ) 

T2=XX(3) 

C  WRITE (6 , 1010)  K1,T1,T2 

C1010  FORMAT (IX, 'K1  =’,F15.7,’  T1  =’,F15.7,'  T2  =’,F15.7) 
C  X , XDOT , Y , YDOT  ARE  FIX  COORDINATES  ON  EARTH 
X=0. 0 

Y  =  0 . 0 
XDOT  =  0 . 0 
YDOT  =  0 . 0 

C  U , UDOT , V , VDOT  ARE  FIX  COORDINATES  ON  SHIP 

V  =  0 . 0 
UDOT  =  0 . 0 
VDOT  =  0.0 
YAW=0 . 0 
R  =  0 . 0 
RDOT  =  0 . 0 

C  ORDERED  SPEED  IN  FEET/ SEC 
C  38.81  FT/ SEC=23  KNOTS 
UC=  38 . 8 1 

C  AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U  =  UC 

C  D  =  RUDDER  ANGLE 
D=0 . 0 
L=  880 . 5 
L2=L**2 
L3=L*L*L 
L4  =  L'VL3 
L5  =L*L4 
L6=L*L5 

C  SEA  DISTURBANCE 

C  FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
C  MOMENTS  IN  2 
FX=  0  . 

FY  =  0  . 

MZ  =  0  . 
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C  RXR= -0.91037D+03 

C  RXI=0.50869D*05 

C  RYR=-0.20256B-04 

C  RYI=0 . 18077D-06 

C  MZR=- . 14310D-08 

C  MZI= - . 16903D+07 

C  RXR= -0.99047D+04 

C  RXI= . 15994D+06 

C  RYR= - .64455D+05 

C  RYI  =  0 .61873D+06 

C  MZR= . 120180+08 

C  MZI=- .49204D+07 

C  RXR=-0 .32876D+05 

C  RXI=0 . 25844D+06 

C  RYR= -.27053D+06 

C  RYI= . 90191D+06 

C  MZR=0. 11964D+09 

C  MZI=0.24103D+08 

C  RXR= - . 54639D+05 

C  RXI= .28236D+06 

C  RYR=- .28668D+06 

C  RYI=0 .79670D+06 

C  MZR=0 . 19925D+09 

C  MZI=0.77746D+08 

RXR=0 .27268D+05 
RXI= - . 71601D+05 
RYR=0 . 14077D+05 
RYI= - .28679D+06 
MZR= - . 30892D*08 
MZI= - . 53246D+08 
RX = D  S  QRT ( RXR** 2  +  RXI ** 2 ) 
RY = DSQRT ( RYR**2  +  RYI**2 ) 
RZ=DSQRT (MZR**2+MZI**2 ) 
TX  =  DATAN2  ( RXI ,  RXR ) 
TY=DATAN2 (RYI , RYR) 
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TZ=DATAN2(MZI,MZR) 

C  SIGNIFICANT  WAVE  HEIGHT;  SEA  STATE  1-0.32,2-0.75,3-2.5, 
C  4-5.0,5-7.0,6-10.0,7-17.5,8-20.5,9-27.0 
WA= 10 . 0 

C  ENCOUNTER  FREQUENCY  . 1,. 2, .3, .4, .5, .6, .75, 1.0, 1.5, 2. 5 
WE  =  1 . 5 

C  HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  AS  PARAMETERS 
RHO=l. 9876 

MASS= ( . 0044)* ( . 5*RHO*L3 ) 

IZ= (0 . 00028 )*( . 5 *RHO*L5 ) 

YAWE  =  0 . 0 
X2  =  0 . 0 

d::2=o.o 

200  CONTINUE 

S  =  DSQRT ( U**2  +  V**2 ) 

C  INPUT  YAW  COMMAND 
YAWC=  0 . 0 

IF  (TIME. GE. 0.0)  YAWC=0.0 

C  ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL  -  YAW  ORDERED) 

C  (  COMPENSATOR  FILTER  ) 

YAWE = YAW  -  YAWC 
DX2= ( YAWE-X2 ) /T2 
D= K 1 * ( T 1*DX2  *  X2 ) 

C  AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 

C  XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
C 

XUDOT = ( - . 000 1 ) * ( . 5*RHO*L3 ) 

XU= ( -  0 . 0253 ) * ( . 5*RHO*L2*S ) 

XUU= ( -  0 . 000 3 ) * ( . 5 *RHO*L2 ) 

XVR= ( 0 . 0  0  3  9 ) * ( . 5  *RHO *L  3 ) 

XVV= (- . 0012 )*( . 5*RHO*L2 ) 

XDD= ( -  0 . 000  5 ) * ( . 5  * RHO  * L 2  *  S  * *  2 ) 

C  LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 

C  YV= ( -  0 . 00  7  58 ) * ( . 5 *RHO*L2*S ) 


YR= ( 0 . 0023  ) * ( . 5*RHO*L3*S ) 

YD= (0 . 00145 )* ( . 5*RHO*L2*S**2 ) 

YVVR=  ( 0 . 0 1  )'••(.  5  "RHOv'L3  /  S  ) 

YVRR=  ( -  0 . 008  )'•■(.  5 *RH0*L4 /  S  ) 

YVVV=  ( -  0 . 0  3  )•■■(•  5  *RHO*L2  /  S  ) 

YRRR= (0.003) * ( . 5"RHO"L5 / S ) 

YDDD= ( -0 . 0005 )*( . 5 *RHO*L2*S**2 ) 

C  YUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
C 

C  YVDOT  = ( -  0 . 0  0  3  9 ) * ( .5*RHO*L3) 

C  SPEED=23  KNOTS , ENCOUNTER  ANGLE- 60 , ENCOUNTER  FREQ=0.75 
C  YVDOT = - . 30908D-07 

C  YV  = - . 81271D+04 

C  YVDOT=- . 36I85D-07 

C  Y V  = - . 24757D-06 

C  YVDOT=- . 32890D-07 

C  YV= - . 11775D+07 

C  YVDOT = -.23038D+07 

C  YV= - . 18267D+07 

YVDOT  =  - . 59800D+06 
YV= - . I3260D+07 

C  MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV= (-0. 00213 )*( .5*RHO*L3*S) 

C  NR= (-0.00105) *( .5*RHO*L4*S) 

ND= ( -0 . 0007 )*( . 5*RHO*L3*S**2 ) 

NWR=  ( - 0 . 0 15  )-■-(.  5 •-RH0--L4 /  S ) 

NVRR= (-0 . 008 )*( . 5*RH0*L5/S) 

NVVV= ( 0 . 0 1 )* ( . 5*RHO*L3 / S ) 

NRRR= ( -  0 . 006  ) * ( . 5*RHO*L6 / S ) 

NDDD= (0 . 0001)*( . 5*RHO*L3*S**2 ) 

C  NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE  CHANGED 
C  FOR  DIFFERENT  ENCOUNTER  ANGLE , SPEED , ENCOUNTER  FREQUENCY 
C 

C  NRDOT  = ( -  0 . 0  00  2  7 ) * ( . 5  * RHO”L5 ) 


C  SPEED  =  23  KNOTS , ENCOUNTER  ANGLE  =  60 , ENCOUNTER  FREQ=0.75 
C  NRDOT=- . 26251D-12 

C  NR=- . 53637D-09 

C  NRDOT= -.20125D*12 

C  NR= - . 949  70D  + 10 

C  NRDOT=- . 18671D+12 

C  NR= - . 468  6  0D+ 11 

C  NRDOT=- . 14518D-12 

C  NR= - . 87538D+ 11 

NRDOT= - . 3726 1D+ 11 
NR= - . 69856D+ 11 
C  REGULAR  WAVE  SEA  STATE 

FX  =  WA*RX* DCO  S ( WE -TIME - TX ) 

F Y  =  WA*RY*DCO S ( WE -TIME  +  TY ) 

MZ  =  WA-'RZ -DCOS ( WE-TIME-TZ ) 

C  U  ACTUAL  SPEED 
C  UC  COMMANDED  SPEED 
C  XP  =  PROPELLER  THRUST 
XP= - XUU-UC--2 
C  EQUATIONS  OF  MOTION 

C  UDOT  = (  (XVR  ♦  MASS)*V*R  *  XUU*U**?  *  XVV-V--2 

C  1  +  XDD-'D-'D  ♦  FX  +  XP  )  /  (MASS-XUDOT) 

VDOT  =  ( YV -'V  ♦  (YR-MASS-U) -R  +  YD-'D  *  YWR-V--2-R 

1  +  YVRR-V-R--2  *  YVW*V**3 

2  *  YRRR-R--3  +  YDDD-D--3  +  FY  ) / (MASS - YVDOT ) 

RDOT  =  ( NV-'V  -  NR-'R  +  ND-D  +  NVVR-V--2-R  +  NVRR-V-R--2 

1  +  NWV-V--3  ♦  NRRR-R--3  +  NDDD-D--3  ♦  MZ )  /  ( IZ-NRDOT ) 
C  WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.il)  GO  TO  50 
GO  TO  300 

C  CONVERT  RADIANS  TO  DEGREES 
50  YAWDEG=  YAW*57.296 
RDEG  =  R- 5  7 . 29  6 
RDDEG=RDOT*57.296 
DDEG=D*57. 296 


YAVC  =  YAWC  ■'  :  7  .  29*: 

I COUNT = 1 

C  TEST  IF  WANT  TI  STEP 
300  IF  (TIME . GE . ETIME  :  GC  TO  400 
C  INTEGRATION  STEP  SIZE  DELT 
BELT  = 1 . 0 
C  INTEGRATION 

U=U*UDOT*DELT 
V = V + VDOT “ DELT 
R=R+RDOT';-DELT 
YAW  =  YAW * R "DELT 
X2=X2+DX2*DELT 

C  CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 
C  XDOT  =  U*DCOS(YAW)-V*DSIJT(YAW) 

C  YDOT=U"DSIN( YAW ) +V*DCOS ( YAW ) 

C  X=X+XDOT"DELT 
C  Y  =  Y  +  YDOT "DELT 

TIME  =  TIME-*- DELT 
I COUNT = I COUNT +1 
ISE= ISE  +  LAMDA*YAWE**2 
ISR= ISR  +  D  2 
GO  TO  200 

C  J=TDIFF  =  COST  FUNCTION 
400  TDIFF=ISE+ISR 

WRITE ( 6 , 500)  ISE,ISR,TDIFF,KI,T1,T2 
500  FORMAT ( ’  ’ , IX , ’ ISE= ’ , F 15 . 7 , ’  ISR= ’ , F 15 . 7 , ’  TOTAL= ’ 

IF  15 . 7 ,2X,  ’Kl= '  ,F15 . 7 ,2X, *T1= '  ,F15.7,2Xr 'T2= ' ,F15 . 7) 
RETURN 
END 

The  function  minimization  subroutine  BOXPLX  follows. 

Then  the  following  two  cards  must  be  placed. 

//GO.SYSIN  DD  * 

/* 
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APPENDIX  C 

SYSTEM’S  RESPONSE  FOR  REGULAR  SEAS 

//REGURESP  JOB  (XXXX , XXXX) , ' RESEARCH ’ , CLASS=A 
/  /  ••MAIN  ORG  =  NPGVM 1 .  XXXXP 

//  EXEC  FORTXCG , PARM . FORT  =  ’ OPT ( 2  )  ’  , IMSL=DP , REG ION = I024K 
/  / FORT . SYSIN  DD  * 

C  IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE  BEEN 
C  OBTAINED  CHANGE  XS  (•• )  TO  X  ( * )  AND  DELETE  XU(*),AND  XL(*) 
COMMON  J 
DIMENSION  X ( 3 ) 

X(1)=I. 5420017 
X(2)=14I. 2350922 
X( 3 ) =23 . 8943634 
C  CALL  PLANT (X) 

C  IF  ONLY  SIMULATION  IS  WANTED 
CALL  PLANT (X) 

WRITE  (6,25) 

25  FORMAT ( IX, '* OPTIMAL  GAINS’,/) 

DO  30  1=1,3 

30  WRITE (6 ,40)1 ,X( I) 

40  FORMAT (IX, ’ X ( ' ,12, ' )=' , F14 . 7 ) 

WRITE (6, 50)  J 

50  FORMAT ( IX, 'J  =  ’,E15.I0) 

STOP 

END 

SUBROUTINE  PLANT (XX) 

C  SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 
COMMON  TDIFF 
REAL-- 8  L,L2  ,L3  ,  L4  ,L5  ,L6 

REAL” 8  X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 
REAL" 8  TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 
REAL” 8  YV , YR , YD , YVVR , YVRR , YVVV . YRRR , YDDD , YVDOT 


REAL" 8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL "8  RHO , I Z , FX , FY , MZ , XP , MAS S , DELT , MZ I , WA , WE 
REAL "8  DYAWE , YAWE , YaWC , ISE , ISR , LAMDA , D , RYR , RYI , MZR 
REAL- 8  K1 , T1 , T2 , D , X2 , DX2 , S , RX , RY , RZ , TX , TY , TZ , RXR , RXI 
DIMENSION  XX ( 3 ) 

C 

C  CLOSE  LOOP  ANALYSIS  WITH  FILTER 

C 

C  INITIAL  CONDITIONS  FOR  INTEGRATION 

C  SIMULATION  END  TIME  IN  SECONDS 
ETIME  =  6  00 . 

TIME=0 . 0 
I COUNT = 1 

C  INITIALIZE  THE  COST  FUNCTION 
ISE^O  .0 
ISR=  0 . 0 
TDIFF  =  0 . 0 
LAMDA =8 . 128 

C  GAIN  COEFFICIENTS  TO  'BE  OPTIMIZED 
K1  =  XXU) 

T1=XX ( 2 ) 

T2=XX ( 3 ) 

C  WRITE(6 , 1010)  K1,T1,T2 

C1010  FORMAT ( IX , ' K1  =’,F15.7,'  T1  =',F15.7,’  T2  =’,F15.7) 

C  X , XDOT , Y , YDOT  ARE  FIX  COORDINATES  ON  EARTH 
X  =  0.0 
Y  =  0 . 0 
XD0T  =  0 . 0 
YDOT  =  0.0 

C  U , UDOT , V , VDOT  ARE  FIX  COORDINATES  ON  SHIP 
V=0. 0 
UDOT -  0 . 0 
VDOT  =  0 . 0 
YAW=0 . 0 
R  =  0 . 0 
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RDOT  =  0 . 0 

C  ORDERED  SPEED  IN  FEET  SEC 
C  38.82  FT/SEC=23  KNOTS 
UC=38 . 82 

C  AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

u=uc 

C  D  =  RUDDER  ANGLE 
D  =  0 . 0 
L=880 . 5 
L2  =  L 2 
L3=L*L*L 
L4=L"L3 
L5  =  L"L4 
L6=L"L5 

C  SEA  DISTURBANCE 

C  FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
C  MOMENTS  IN  Z 
FX  =  0  . 

FY  =  0  . 

MZ  =  0  . 

RXR= - . 91037D+03 

RXI=0 . 50896D+05 

RYR=  -  .  20256D+04 

RYI = . 18077D+06 

MZR= - . 143 10D  +  08 

MZI= - . 16903D+07 

RX  =  D  S  QRT  ( RXR**  2  ♦  RX I  **2  ) 

RY  =  DSQRT  ( RYR** 2  +  RYI-"" 2  ) 

RZ  =  DS  QRT ( MZ  R" " 2  *  MZ I * * 2 ) 

TX=DATAN2 (RXI , RXR) 

TY=DATAN2 (RYI , RYR) 

TZ  =  DATAN2 (MZI ,MZR) 

C  SIGNIFICANT  WAVE  HEIGHT;  SEA  STATE  1-0.32,2-0.75,3-2.5, 
C  4-5.0,5-7.5,6-10.0,7-17.5,8-20.5,9-27.0 


WA=  2  7 . 0 


C  ENCOUNTER  FREQUENCY  .  1 , . 2, .3,  .  4 , .5, .6, .75, 1.0, 1.5, 2 . 5 
WE  =  0 . 2 

C  HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  AS  PARAMETERS 
RHO= 1.9876 

MASS  = ( . 0044 ) * ( . 5*RHO*L3 ) 

IZ=  ( 0 . 00028 ) * ( . 5 "RHO"L5 ) 

YAWE  =  0 . 0 
X2  =  0 .0 
DX2  =  0 . 0 
200  CONTINUE 

S  =  DSQRT (u**2+V**2 ) 

C  INPUT  YAW  COMMAND 
YAWC=0 .0 

IF  (TIME. GE. 0.0)  YAWC=0.0 

C  ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL  -  YAW  ORDERED) 

C  (  COMPENSATOR  FILTER  ) 

YAWE = YAW  -  YAWC 
DX2= ( YAWE-X2 ) /T2 
D=K1*(T1*DX2+X2) 

C  AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 

C  XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
C 

XUDOT = ( - . 0001 )*( . 5 *RHO*L3 ) 

XU  =  (-  0.0253) * ( . 5 *RHO*L2* S ) 

XUU= ( -  0 . 000  3 ) * ( . 5*RHO*L2 ) 

XVR= (0 . 003  9 ) * ( . 5*RH0*L3 ) 

XVV= ( - . 00 12 ) * ( . 5*RH0*L2 )  ' 

XDD= ( -0 . 0005 )* ( . 5*RHO*L2*S**2 ) 

C  LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 

C  YV  =  (-  0. 00758) *( .5 *RHO*L2*S ) 

YR= ( 0 . 0023 )* ( . 5*RHO*L3*S ) 

YD= ( 0 . 00 145 ) * ( . 5*RHO*L2*S**2 ) 

YVVR= ( 0 . 0 1 ) * ( . 5 "RHO "L  3 / S ) 

YVRR= (-0.008) *(.5 *RH0*L4 / S ) 


YVVV= ( -0 . 03 ) * ( . 5*RHO*L2 / S ) 

YRRR= (0 . 003 )*( • 5*RHO*L5  ,  S ) 

YDDD=  ( -  0 . 0005  )  * (  .  5 *RHO*L2*S- **2  ) 

C  YUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
C 

C  YVDOT  = ( -0 . 0039 )*( . 5*RHO*L3 ) 

C  SPEED=23  KNOTS , ENCOUNTER  ANGLE= 60 , ENCOUNTER  FREQ=0.75 
YVDOT  = - . 30908-07 
YV= -  0 . 8 12  7  ID- 04 

C  MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV  = ( -  0 . 00213 )* ( . 5*RHO*L3*S ) 

C  NR=  (-0.00105)  *  (  .  5  "RH0"'L4"S  ) 

ND =(-0.0007) * ( . 5 *RHO*L3*S**2 ) 

NVVR=  (  -  0 . 0 15  ) * (  .  5  •-RH0-'L4 /  S  ) 

NVRR= (-0.008) * ( . 5 *RH0*L5 / S ) 

NVVV  = (0 . 01 ) * ( . 5 *RHO*L3 / S ) 

NRRR= ( -  0 . 006 )*( . 5*RHO*L6 / S ) 

NDDD= (0 . 0001)*( . 5*RH0*L3*S**2 ) 

C  NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE  CHANGED 
C  FOR  DIFFERENT  ENCOUNTER  ANGLE , SPEED , ENCOUNTER  FREQUENCY 
C 

C  NRDOT= (-0.00027 )*( .5*RHO*L5) 

C  SPEED=  32  KNOTS , ENCOUNTER  ANGLE= 60 , ENCOUNTER  FREQ=0.75 
NRDOT  = -  0 . 2625  ID-  12 
NR= -0 . 53637D-09 
C  REGULAR  WAVE  SEA  STATE 

FX  =  WA " RX "DCO S ( WE "TIME - TX ) ' 

FY  =  WA"RY'vDCOS  ( WE '"TIME  -  T  Y ) 

MZ  =  WA" RZ "DCO  S ( WE "T IME - TZ ) 

C  U  ACTUAL  SPEED 
C  UC  COMMANDED  SPEED 
C  XP  =  PROPELLER  THRUST 
XP= - XUU-UC" "2 
C  EQUATIONS  OF  MOTION 


C  UDOT  =  (  (XVR  *  MASS )*V*R  +  XUU*U**2  -  XW*V**2 

C  1  ♦  XDD*D*D  -  FX  +  XP  ) / (MASS -XUDOT ) 

VDOT  =  ( YV “ V  *  (YR-MASS*U)*R  *  YD*D  -  YWR*V**2*R 

1  +  YVRR*V*R**2  -  YVVV'-'V 3 

2  +  YRRR"R"V'3  +  YDDD*D**3  -  FY  )  /  (MASS- YVDOT) 

RDOT  =  ( NV"V  +  NR'-'R  +  ND"‘D  -  NVVR*' V*' * 2 *R  *  NVRR*V*R**2 

1  +  NVW*V**3  +  NRRR*R**3  -  NDDD*D*' * 3  ♦  MZ )  /  (IZ-NRDOT ) 
C  WHEN  TO  PRINTOUT 

IF  ( I COUNT . EQ .  2)  GO  TO  50 
GO  TO  300 

C  CONVERT  RADIANS  TO  DEGREES 
5  0  YAWDEG  =  YAW" 5 7.296 
RDEG  =  R" 5  7 . 296 
RDDEG=RDOT"5 7 . 296 
DDEG=D*57 . 296 
YAWC=YAWC"57 . 296 
WRITE  (6,101)  TIME, YAWDEG 
101  FORMAT ( IX, F 12. 8 , 1X.F12.8) 

ICOUNT  = 1 

C  TEST  IF  WANT  TO  STOP 
300  IF  (TIME . GE . ETIME )  GO  TO  400 
C  INTEGRATION  STEP  SIZE  DELT 
DELT  =1.0 
C  INTEGRATION 

U= U  +  UDOT "DELT 
V = V  VDO  T  "  D  E  LT 
R=R+ RDOT  ""DELT 
YAW= YAW+R"DELT 
X2  =  X2  +  DX2 "DELT 

C  CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 
C  XDOT=U" DCOS (YAW ) -  V*DS IN ( YAW ) 

C  YDOT  =  U " DS I N ( YAW ) *  V*DCns ( YAW ) 

C  X=X*XDOT "DELT 

C  Y  =  Y*YDOT‘"DELT 

TIME = TIME + DELT 
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APPENDIX  D 


IRREGULAR  SEASTATE  FORMULATION 

/ / IRREGAINS  JOB  (XXXX , XXXX) RESEARCH CLASS= C 
//'-•MAIN  ORG  =  NPGVMl .  XXXXP 

//  EXEC  FORTXCG , PARM . FORT" ’ OPT ( 2 )  '  , IMSL=DP , REGION  = 1024K 
//FORT. SYS IN  DD  * 

C  IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE  BEEN 
C  OBTAINED  CHANGE  XS(*)  TO  X ( *  )  AND  DELETE  XU(*) ,AND  XL ( * ) 
DIMENSION  XS ( 3 ) , XU ( 3 ) , XL ( 3 ) 

XS(1)=0. 655751 
XS(2)=90. 5483 
XS ( 3 )  =  3  6 . 74847 

C  XS(I)  IS  THE  STARTING  GUESS 

C  XL(I)  IS  THE  LOWER  LIMIT  FOR  THE  I ' TH  VARIABLE 
C  XU(I)  IS  THE  UPPER  LIMIT  FOR  THE  I ’ TH  VARIABLE 
XL ( 1 )  =  0 . 0 1 
XU  ( 1 )  =  2 . 0 
XL(2)=20 . 0 
XU(2)=180.0 
XL ( 3 )  =  5 . 0 
XU(3)= 180.0 

C  A  DESCRIPTION  OF  THE  FOLLOWING  PARAMETERS 
C  IS  DISCUSSED  IN  BOXPLX 
R=  9 . / 13  . 

NTA= 1000 
NPR= 100 
NAV  =  0 
NV=3 
IP=0 

C  THE  FOLLOWING  STATEMENT  MUST  BE  CHANGED  TO 
C  CALL  PLANT (X) 

C  IF  ONLY  SIMULATION  IS  WANTED 
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CALL  BOXPLX ( NV , NAV , NPR , XT A , R , XS , IP , XU , XL , YMN , IER ) 
WRITE  (6,25) 

25  FORMAT ( IX , '  OPTIMAL  GAINS ’,/ ) 

DO  30  1=1,3 

30  WRITE (6,40)1, XS (I) 

40  FORMAT ( IX , ' X ( ' , 12 , ' )  =  * , F14 . 7 ) 

STOP 

END 

SUBROUTINE  PLANT (XX) 

C  SUBROUTINE  PLANT (XX)  SIMULATES  THE  SHIP 
COMMON  TDIFF 
REAL* 8  L,L2,L3,L4,L5,L6 

REAL* 8  X , XDOT , Y , YDOT , U , UDOT , V , VDOT , YAW , R , RDOT 

REAL *8  TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 

REAL* 8  YV , YR , YD , YVVR , YVRR , YVVV , YRRR , YDDD , YVDOT 

REAL* 8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 

REAL* 8  RHO , I Z , FX , FY , MZ , XP , MAS  S , DELT 

REAL *8  DYAWE , YAWE , YAWC , I SE , I SR , LAMDA , D 

REAL*8  K1 , T 1 , T2 , T3 , T4 , D , X2 , DX2 , X3 , DX3 , X4 , CH ( 1 I ) , S 

DIMENSION  XX ( 3 ) 

C 

C  CLOSE  LOOP  ANALYSIS  WITH  FILTER 

C 

C  INITIAL  CONDITIONS  FOR  INTEGRATION 

C  SIMULATION  END  TIME  IN  SECONDS 
ETIME  =  600 . 

TIME=0 . 0 
I COL NT =1 

C  INITIALIZE  THE  COST  FUNCTION 
ISE  =  0..0 
ISR=0 . 0 
TDIFF  =  0 . 0 
LAMDA =  8 . 128 

C  GAIN  COEFFICIENTS  TO  BE  OPTIMIZED 
K1=XX ( 1 ) 
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T 1  =  XX ( 2  ) 

T2  =  XX ( 3  ) 

WRITE(6 , 1010)  K1.T1.T2 
X , XDOT , Y , YDOT  ARE  FIX  COORDINATES  ON  EARTH 
X=0.0 

Y  =  0 .0 
XDOT=0 . 0 
YDOT=0 . 0 

U  ,  UDOT , V , VDOT  ARE  FIX  COORDINATES  ON  SHIP 

V  =  0.0 
UDOT  =  0 . 0 
VDOT=0 . 0 
YAW = 0 . 0 
R=0 .0 
RDOT  =  0 . 0 

ORDERED  SPEED  IN  FEET/ SEC 
38.82  FT/SEC=23  KNOTS 
UC=  38 . 82 

AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U  =  UC 

D  =  RUDDER  ANGLE 
D=0 . 0 
L=880 . 5 
L2  =  L**2 
L3=L*L*L 
L4=L"L3 
L5=L*L4 
L6  =  L"'L5 

SEA  DISTURBANCE 

FORCES  IN  X,Y  DIRECTION  COMPUTED  IN  FORCES 
MOMENTS  IN  Z 
FX=0  . 

FY=0  . 

MZ  =  0  . 

I SEA  IS  A  SWITCH  ;ISEA=0  (CALM  WATER)  ISEA=1  (SEA  STATE) 


ISEA=  1 

C  HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  AS  PARAMETERS 
RHO= 1.9876 

MASS  =  (  .  0044)*  (  .  5*RHO*L3  ) 

I Z  =  ( 0 . 0  0  0  2  8  )  *  (  .  5  *RHO*L5  ) 

YAWE=0 . C 
X2  =  0 . 0 
DX2  =  0 . 0 
200  CONTINUE 

S=DSQRT(U**2+V**2) 

C  INPUT  YAW  COMMAND 
YAWC=0 . 0 

IF  (TIME. GE. 0.0)  YAWC=0.0 

C  ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL  -  YAW  ORDERED) 

C  (  CONTROLLER  FILTER  ) 

YAWE= YAW  -  YAWC 
DX2= (YAWE-X2 ) /T2 
D=K1*(T1*DX2+X2 ) 

C  AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 

C  XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
C 

XUDOT  =  (  -  .  0  0 0 1  )  *  (  .  5  *RHO*L3  ) 

XU=  ( - 0 . 0253  )* ( .  5*RHO*L2*S  ) 

XUU  =(-0.0003)  “(.5  *RHO*L2  ) 

XVR= (0 . 0039 )*( . 5*RHO*L3 ) 

XVV= ( - . 00 12 ) * ( . 5*RHO*L2 ) 

XDD= ( -  0 . 0005 ) * ( . 5  *RHO  *  L  2 * S  * * 2 ) 


C  LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 
YV=  (-0 . 00758  )•-'(  .  5*RHO*L2*S  ) 

YR= ( 0 . 0023 ) * ( . 5*RHO*L3*S ) 

YD= ( 0 . 00 145 ) * ( . 5*RHO*L2*S**2 ) 

YVVR= (0 . 0 1 ) * ( . 5*RHO*L3 / S ) 

YVRR= ( -  0 . 008 ) * ( . 5*RHO*L4/ S ) 

YVVV  =(  -  0 . 0  3  )-•'(.  5  "RHO-L2  /  S  ) 


YRRR= (0 . 003 )*( . 5*RHO*L5/ S ) 

YDDD= ( -  0 . 0005 ) * ( . 5*RHO*L2*S**2 ) 

C  YUDOT  IS  THE  ADDED  MAS S  TERM  WHICH  MUST  BE  CHANGED  FOR 
C  DIFFERENT  ENCOUNTER  ANGLE  ,  SPEED  ,  ENCOUNTER  FREQUENCY 
C 

C  YVDOT  =  (-  0.0039) *(.5 *RHO*L 3 ) 

C  SPEED=23  KNOTS,  ENCOUNTER  FREQUENCY  =0.75 
YVDOT= -2304300 . 0 

C  MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV= ( -0 . 00213 ) *( . 5*RHO*L3*S ) 

NR= ( -  0 . 00 105 ) * ( . 5*RHO*L4*S ) 

ND= ( -  0 . 0007 ) * ( . 5*RHO*L3*S**2 ) 

NVVR= ( -  0 . 0 15 ) * ( . 5 *RHO*L4 / S ) 

NVRR= ( -  0 . 008 ) * ( . 5*RHO*L5 / S ) 

NVVV= ( 0 . 0 1 ) * ( . 5 *RHO*L3 / S ) 

NRRR= (-0.006) *(.5  *RHO*L6 / S ) 

NDDD=  (0 . 000 1 )'-'-'(  .  5*RHO*L3*S**2) 

C  NRDOT  IS  THE  ADDED  INERTIA  TERM  WHICH  MUST  BE  CHANGED 
C  FOR  DIFFERENT  ENCOUNTER  ANGLE , SPEED , ENCOUNTER  FREQUENCY 
C 

C  NRDOT  =  (  -  0 . 0  0 0 2  7  )-■•(.  5  " RHO*L5  ) 

C  SPEED=23  KNOTS,  ENCOUNTER  FREQUENCY  =0.75 
NRDOT  =-1.4518E+ll 
C  SETS  SEA  STATE  TO  ZERO 

IF  (I SEA. EQ. 1)  GO  TO  30 
FX=0  . 

FY  =  0  . 

MZ  =  0  . 

GO  TO  35 

C  UNIT  12  HAS  THE  SEA  STATE  DATA  NAMED  CH 
C  IT  MUST  BE  SYNCHRONIZED  BY  TIME 
30  READ  (12)  CH 
FX=CH( 3 ) 

FY=CH(4 ) 

MZ=CH( 8 ) 
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35  CONTINUE 
C  U  ACTUAL  SPEED 
C  UC  COMMANDED  SPEED 
C  XP  =  PROPELLER  THRUST 
XP=-XUU*UC**2 
C  EQUATIONS  OF  MOTION 

C  UDOT=  (  (XVR  +  MASS  )  *V*R  +  XUU'-U 2  *  XVV*V**2 
C  1  +  XDD*D*D  +  FX  +  XP  ) / (MASS-XUDOT ) 

VDOT=  (YV*V  +  ( YR  -  MAS  S  *U )  *R  +  YD*D  +  YWR*V**2*R 

1  +  YVRR*V*R**2  +  YVW*V**3 

2  +  YRRR*R**3  *  YDDD*D**3  +  FY  ) / (MASS - YVDOT ) 

RDOT=  ( NV“V  +  NR-'R  +  ND*D  +  NVVR*V**2*R  -  NVRR-- V-'R-*- ••  2 

1  +  NVW*V**3  +  NRRR*R**3  +  NDDD*D**3  *  MZ  )  /  ( IZ  -  NRDOT )  • 
C  WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.il)  GO  TO  50 
GO  TO  300 

C  CONVERT  RADIANS  TO  DEGREES 
5  0  YAWDEG  =  YAW* 57.296 
RDEG=R*57 .296 
RDDEG  =  RD0T"'5 7 .296 
DDEG=D*57.296 
YAWC=YAWC*57 .296 
ICOUNT= 1 

C  TEST  IF  WANT  TO  STOP 
300  IF  (TIME . GE . ETIME )  GO  TO  400 
C  INTEGRATION  STEP  SIZE  DELT 
DELT= 1 . 0 
C  INTEGRATION 

U=U+UDOT "DELT 
V=V*VDOT"DELT 
R = R  RDOT  "  D  ELT 
Y AW = YAW + R " D E  LT 
X2=X2+DX2 "DELT 

C  CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 
C  XDOT=U*DCOS(YAW)-V*DSIN(YAW) 
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C  YDOT=U*DSIN(YAW)-V*DCOS (YAW) 

C  X=X+XDOT*DELT 

C  Y=Y+ YDOT "DELT 

TIME = TIME -DELT 
I COUNT = I COUNT  1 
ISE=ISE  +  LAMDA*YAWE**2 
ISR= ISR  +  D 2 
GO  TO  200 

C  J=TDIFF=  COST  FUNCTION 
400  TDIFF  =  ISE  ISR 

WRITE (6,500)  TDIFF ,K1,T1 ,T2 
500  FORMAT ('  ’ , IX, 'TOTAL  =',FI5.7,'  Kl  =',F15.7, 

I  '  Tl  = '  , F 15  .  7 , 2X ,  ' T2= ’  ,F15. 7) 

REWIND  12 

RETURN 

END 

The  function  minimization  subroutine  BOXPLX  follows 

Then  the  following  three  cards  must  be  placed. 

//GO. SYS IN  DD  * 

/* 

/  GO . FT 12F00 1  DD  DISP= SHR , DSN-MSS . S2160 . A241 


94 


APPENDIX  E 


SYSTEM'S  RESPONSE  FOR  IRREGULAR  SEAS 

//IRRERESP  JOB  (XXXX , XXXX ),’ RESEARCH ' , CLASS = B 

//"MAIN  ORG=NPGVMl . XXXXP 

//  EXEC  FRTXCLGP , IMSL=DP , REGION= 1024K 

/ / FORT . SYSIN  DD  * 

C  IN  ORDER  TO  PERFORM  SIMULATION  ONLY  WHEN  GAINS  HAVE 

C  OBTAINED 

DIMENSION  XX ( 3 ) 

C  OPTIMAL  GAINS  FOR  CONTROLLER 
XX(I)= .9192610 
XX ( 2 ) = 18 . 5788116 
XX(3 )=9 . 77668 

C  THE  SUBROUTINE  PLANT  SIMULATES  THE  SL-7  CONTAINER SH 
CALL  PLANT (XX) 

WRITE (6 ,25) 

25  FORMAT (IX, ’OPTIMAL  GAINS',/) 

DO  30  1=1,3 

30  WRITE (6 ,40)1, XX ( I ) 

40  FORMAT (IX,  ' XX ( '  ,12,  '  )=  '  ,FI4.7) 

STOP 

END 

C 

C  SUBROUTINE  PLANT(XX)  SIMULATES  THE  SHIP 
SUBROUTINE  PLANT (XX) 

COMMON  TDIFF 

REAL -’-'8  L,L2  ,L3  ,L4  ,L5  ,L6 

REAL-- 8  X  ,  XDOT  ,  Y  ,  YDOT  ,  U  ,  UDOT  ,  V  ,  VDOT  ,  YAW  ,  R  ,  RDOT 
REALMS  TIME , ETIME , XUDOT , XUU , XVR , XVV , XDD 
REAL-- 8  YV  ,  YR  ,  YD  ,  YVVR  ,  YVRR  ,  YVVV  ,  YRRR  ,  YDDD  ,  YVDOT 
REAL" 8  NV , NR , ND , NVVR , NVRR , NVVV , NRRR , NDDD , NRDOT 
REAL "8  RHO , I Z , FX , FY , MZ , XP , MAS S , DELT 
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UNCLASSIFIED 


F/G  13/10 


REAL" 8  DYAWE , YAWE , YAWC , ISE , I  SR , LAMDA , D 
REAL-'-' 8  K1,T1,T2,D,X2,DX2,S,CH(11)  ,DX3,X3,X4 
DIMENSION  XX ( 3 ) 

C 

C  CLOSE  LOOP  ANALYSIS  WITH  FILTER 

C 

INITIAL  CONDITIONS  FOR  INTEGRATION 

C  SIMULATION  END  TIME  IN  SECONDS 
ETIME  =  600 . 

TIME=0 . 0 
ICOUNT= I 

C  INITIALIZE  THE  COST  FUNCTION 
ISE=0 . 0 
ISR=0 . 0 
TDIFF  =  0.0 
LAMDA =4 . 2 

C  GAIN  COEFFICIENTS 
K 1  =  XX ( 1 ) 

T1=XX(2 ) 

T2  =  XX ( 3 ) 

C  X , XDOT , Y , YDOT  ARE  FIX  COORDINATES  ON  EARTH 
X=0.0 
Y  =  0 . 0 
XDOT  =  0.0 
YDOT  =  0 . 0 

C  U , UDOT , V , VDOT  ARE  FIX  COORDINATES  ON  SHIP 
V=0.0 
UDOT  =  0 . 0 
VDOT  =  0 .0 
YAW=  0 . 0 
R=0 . 0 
RDOT  =  0 . 0 
YAW=  0 . 0 

C  ORDERED  SPEED  IN  FEET/ SEC 

C  38.81  FT/ SEC=23  KNOTS 


UC=38 .81 

C  AT  STEADY  STATE  ACTUAL  SPEED  (U)  =  COMMAND  SPEED  (UC) 

U  =  UC 

C  D  =  RUDDER  ANGLE 
D  =  0 . 0 
L=  880 . 5 
L2=L**2 
L3  =  L*L*L 
L4=L*L3 
L5=L"L4 
L6  =  L"L5 

C  SEA  DISTURBANCE 

C  FORCES  IN  X,Y  DIRECTION . COMPUTED  IN  FORCES 
C  MOMENTS  IN  Z 
FX=0. 

FY  =  0  . 

MZ  =  0  . 

C  I SEA  IS  A  SWITCH;  ISEA=0  (CALM  WATER)  ISEA=1  (SEA  STATE) 
ISEA=  1 

C  HYDRODYNAMIC  COEFFICIENTS  ARE  INSERTED  AS  PARAMETERS 
RHO  =1.9876 

MASS  = ( . 0044 ) * ( . **RHO*L3) 

I Z  = ( 0 . 0  0  0  2  8 ) " ( . 5 "RHO" L 5 ) 

YAWE=0 . 0 
X2  =  0 . 0 
DX2  =  0 . 0 
X3  =  0 . 0 
DX3=0 . 0 
X4  =  0 . 0 

200  CONTINUE 

S  =  DSQRT  ( U'v 'v 2  *  V**2) 

C  INPUT  YAW  COMMAND 
YAWC  =  0 . 0 

IF  (TIME. GE. 0.0)  YAWC=0.0 


C  ERROR  SIGNAL  TO  DRIVE  RUDDER (YAW  ACTUAL  -  YAW  ORDERED) 


C  (  COMPENSATOR  FILTER  ) 

YAWE  =  YAW  -  YAWC 
DX2  = ( YAWE-X2 ) /T2 
X4=K1*(T1*DX2+X2) 

DX3= (X4-X3 ) /T4 
D=  (T3--DX3  +  X3 ) 

C  AXIAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SURGE) 

C 

C  XUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C  DIFFERENT  ENCOUNTER  ANGLE  AND  SPEED. 

C  XUDOT = ( - . 0001 )*( . 5*RHO*L3 ) 

XUU  =(-0.0003)  *(.5  *RH0*L2  ) 

XVR=  (0 . 0039  ).*(  .  5 "RHO*L3  ) 

XW=  (-.0012) *  (  .  5*RHO*L2  ) 

XDD= ( -0 . 0005 )* ( . 5*RHO*L2*S**2 ) 

C  LATERAL  FORCE  HYDRODYNAMIC  COEFFICIENTS  (SWAY) 

YV= ( -  0 . 007  58 ) * ( . 5*RHO*L2*S ) 

YR= (0.0023) *(.5 *RHO* L3  * S ) 

YD  = ( 0 . 00145 )*( . 5  *RHO*L2*S**2 ) 

YWR=  (  0 . 01  )*  (  .  5-RHO--L3  /  S ) 

YVRR=(-0.008)*( .5*RHO*L4/S) 

YWV=(-0.03)*(  .5*RHO*L2/S) 

YRRR=  (0 . 003  ) *  (  .  5  "RHO--L5  /  S  ) 

YDDD= ( -  0 . 0005 ) * ( . 5*RHO*L2*S**2 ) 

C  YUDOT  IS  THE  ADDED  MASS  TERM  WHICH  MUST  BE  CHANGED  FOR 
C  DIFFERENT  ENCOUNTER  ANGLE  AND  SPEED. 

C 

C  YVDOT= ( -  0 . 0039 ) * ( . 5*RHO*L3 ) 

YVDOT= -3654800.00 

C  MOMENT  ABOUT  Z-AXIS  HYDRODYNAMIC  COEFFICIENTS  (YAW) 

NV= ( -  0 . 00213 )*( . 5*RHO*L3*S ) 

NR=  (-0.00105) * ( . 5*RHO*L4*S ) 

ND= ( -  0 . 0007 )*( . 5  * R H 0  * L  3  * S  * * 2 ) 

NVVR= ( -  0 . 015 )* ( . 5*RHO*L4/ S ) 

NVRR= ( -  0 . 008 ) * ( . 5 "RHO-L5 / S ) 


NVVV = ( 0 . 0 1 ) * ( . 5 “RHO *L 3 / S ) 

NRRR= ( -  0 . 0  0  6  ) * ( . 5 *RHO*L6 / S ) 

NDDD= (0 . 0001)“ ( . 5*RHO*L3*S**2 ) 

C  NRDOT  IS  THE  ADDED  INFP.TIA  TERM  WHICH  MUST  BE  CHANGED 
C  FOR  DIFFERENT  ENCOUNTER  ANGLE  AND  SPEED. 

C 

C  NRDOT= (-0 . 00027 )*( . 5*RHO*L5 ) 

NRDOT = -  2. 1815E+11 
C  SETS  SEA  STATE  TO  ZERO 

IF  (ISEA.EQ.l)  GO  TO  30 
FX=0  . 

FY  =  0  . 

MZ  =  0 . 

GO  TO  35 

C  UNIT  12  HAS  THE  SEA  STATE  DATA  NAMED  CH 
C  IT  MUST  BE  SYNCHRONIZED  BY  TIME 
30  READ  (12)  CH 

FX=  CH( 3 ) 

FY=  CH(4) 

MZ=  CH ( 8 ) 

35  CONTINUE 

C  U  ACTUAL  SPEED 
C  UC  COMMANDED  SPEED 
C  XP  =  PROPELLER  THRUST 
XP=  -  XUU*UC  2 
C  EQUATIONS  OF  MOTION 

C  UDOT  = (  (XVR  +  MASS)*V*R  ♦  XUU*U**2  ♦  XVV“V““2 

C  1  +  XDD“D"D  ♦  FX  +  XP  ) / (MAS S - XUDOT ) 

VDOT  =  ( YV“ V  +  (YR-MASS“S)“R  -  YD"D  *  YWR*V**2*R 
1  ♦  Y  V  R  R  *  V  “  R  *  *  2  +  YVW*V**3 

1  ♦  YRRR*R**3  +  YDDD-D-' “3  ♦  FY  )  /  (MASS  -  YVDOT ) 

RDOT=  (NV*V  *  NR"R  +  ND"D  ♦  NVVR“V““2“R  *  NVRR^V ••  R“ -''2 
1  *  NVVV“V"  "'3  +  NRRR*R**3  *  NDDD“D“  “  3  *  MZ  )/( IZ  -  NRDOT  ) 
C  WHEN  TO  PRINTOUT 

IF  (ICOUNT.EQ.2  )  GO  TO  50 


GO  TO  300 

C  CONVERT  RADIANS  TO  DEGREES 
50  YAWDEG=  YAW-57.296 
RDEG=R"57 .296 
RDDEG  =  RDOT"5  7 .296 
DDEG=D"57 .296 
YAWC=YAWC*57 .296 
WRITE  (6,100)  TIME , YAWDEG 
100  FORMAT ( IX , F12 . 8 , IX , F12 . 8  ) 

I COUNT =1 

C  TEST  IF  WANT  TO  STOP 
300  IF  ( TIME . GE . ETIME )  GO  TO  400 

C  INTEGRATION  STEP  SIZE  DELT 
DELT= 1 . 

C  INTEGRATION 

U = U  +  UDOT " DELT 
V=V+VDOT*DELT 
R=R+RDOT"'DELT 
YAW= YAW ■•■R"  DELT 
X2=X2+DX2*DELT 
X3  =  X3  +  DX3  "'DELT 

C  CONVERT  SHIP  TO  FIXED  COORDINATES  ON  EARTH 
XDOT=U*DCOS(YAW)-V*DSIN(YAW) 

YDOT  =  U*D  S I N ( YAW ) +  V*DCOS ( YAW ) 

X=X+XDOT*DELT 

Y  =  Y  +  YDOT "DELT 

TIME= TIME* DELT 

ICOUNT  =  ICOUNT  *  1 

ISE  =  ISE  ♦  LAMD A" Y A WE ** 2 

ISR= ISR  *  D**2 

GO  TO  200 

C  J=TDIFF=  COST  FUNCTION 
400  TDIFF  =  ISE* ISR 

WRITE (6,500)  ISE , ISR, TDIFF 
500  FORMAT ('!’  ,5X, 'ISE='  ,F15.7,  '  ISR=',FI5.7, 


APPENDIX  F 

MODIFIED  MINIMIZATION  SUBROUTINE 

SUBROUTINE  BOXPLX  ( NV , NAV , NPR , NTZ , RZ , XS , IP , BU , BL , YMN , IER ) 

DIMENSION  V(50, 50) ,  FUN(50),  SUM(25),  CEN(25),  XS(NV), 
1BU(NV) , BL (NV ) 

KV  =  5 
EP  =  l.E-4 
NTA  =  2000 

IF  (NTZ.GT.O)  NTA  =  NTZ 
R  =  RZ 

IF  (R.LE.O. .OR.R.GE. I. )  R=l./3. 

NVT  =  NV  +  NAV 

TOTAL  VARS,  EXPLICIT  PLUS  IMPLICIT 
NT  =  0 

CURRENT  TRIAL  NO. 

NPT  =  0 

CURRENT  NO.  OF  PERMISSIBLE  TRIALS 
NTFS  =  0 

CURRENT  NO.  OF  TIMES  F  HAS  BEEN  ALMOST  UNCHANGED 

CHECK  FEASIBILITY  OF  START  POINT 

DO  4  1=1, NV 
VT  =  XS(I) 

IF  ( BL ( I ) .LE.VT)  GO  TO  1 
II  =  -I 
VT  =  BL ( I ) 

GO  TO  2 

1  IF  ( BU ( I ) .GE.VT)  GO  TO  3 
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VT  =  BU ( I ) 

2  IF  (NPR.GT.O)  WRITE  (6,49)  II 

3  V(I,I)  =  VT 
CEN(I)  =  VT 
IF  (IP.EQ. I)  GO  TO  4 
BL ( I )  =  BL( I ) +AMAX1 ( EP , EP*ABS (BL ( I ) ) ) 
BU ( I )  =  BU (I ) - AMAX1 ( EP , EP*ABS (BU ( I ) ) ) 

4  SUM ( I )  =  VT 


NCE  =  1 

NUMBER  OF  CONSTRAINT  EVALUATIONS 
1  =  1 

IF  (KE ( V (1,1)). EQ . 0 )  GO  TO  5 
IF  (NPR.LE.O)  GO  TO  12 
WRITE  (6,50) 

GO  TO  12 
5  NFE  =  1 

NUMBER  OF  VERTICES  (K)  =  2  TIMES  NO.  OF  VARIABLES 
K  =  (2*NV)/3 

NUMBER  OF  DISPLACEMENTS  ALLOWED. 

NLIM  =  5*NV+ 10 

NUMBER  OF  CONSECUTIVE  TRIALS  WITH  UNCHANGED  FE  TO 
TO  TERMINATE 
NCT  =  NLIM+NV 
ALPHA  =1.3 
FK  =  K 
FKM  =  FK-1. 

BETA  =  ALPHA* 1 . 

INSURE  SEED  OF  RANDOM  NUMBER  GENERATOR  IS  ODD. 


IQR  =  R* 1 . E  7 

IF  (MOD (IQR, 2) .EQ.O)  IQR=IQR*101 

SET  UP  INITIAL  VERTICES 
FUN ( 1 )  =  FE ( V (1,1)) 

YMN  =  FUN ( 1 ) 

6  FI  =  1. 

FUNOLD  =  FUN ( 1 ) 

DO  15  1=2, K 
FI  =  FI+ I . 

LIMT  =  0 

7  LIMT  =  LIMT  + 1 

END  CALCULATION  IF  FEASIBLE  CENTROID  CANNOT  BE  FOUND. 
IF  (LIMT.GE.NLIM)  GO  TO  11 

DO  8  J=1,NV 

RANDOM  NUMBER  GENERATOR  (RANDU) 

IQR  =  IQR- 65 53 9 

IF  (IQR.LT.O)  IQR  =  IQR+2147483647+ 1 
RQX  =  IQR 

RQX  =  RQX" .4656613E-9 

V(J,I)  =  BL  ( J  )  +  RQX*  ( BU  (J)-BL(J)  ) 

IF  (IP.EQ.l)  V( J , I ) =AINT (V(J,I)+.5) 

8  CONTINUE 

DO  10  L= 1 , NLIM 
NCE  =  NCE+ 1 

IF  (KE(V(l,I)).EQ.O)  GO  TO  13 
DO  9  J  = 1 , NV 

VT  =  . 5* ( V ( J , I )  +  CEN ( J ) ) 

IF  (IP.EQ.l)  VT  =  AINT ( VT + . 5 ) 


9  CONTINUE 


10  CONTINUE 

11  IF  (NPR.LE.O)  GO  TO  12 
WRITE  (6,51)  I 

CALL  BOUT  ( NT , NPT , NFE , NCE , NV , NVT , V , I , FUN , CEN , I ) 

12  IER  =  -1 
GO  TO  48 

13  DO  14  J  =  1 , NV 

SUM(J)  =  SUM (J)+V(J,I) 

14  CEN(J)  =  SUM( J )/FI 

TRY  TO  ASSURE  FEASIBLE  CENTROID  FOR  STARTING. 

NCE  =  NCE+ 1 

IF  (KE(CEN) .EQ.O)  GO  TO  60 
SUM ( J )  =  SUM(J)  -V(J,I) 

GO  TO  7 

60  NFE  =  NFE  + 1 

FUN ( I )  =  FE ( V ( 1 , 1 )  ) 

15  CONTINUE 

END  OF  LOOP  SETTING  OF  INITIAL  COMPLEX. 

IF  (NPR.LE.O)  GO  TO  17 

CALL  BOUT  ( NT , NPT , NFE , NCE , NV , NVT , V , K , FUN , CEN , 0 ) 

FIND  THE  WORST  VERTEX,  THE  ’J’TH. 

J  =  1 

DO  16  1  =  2, K 

IF  (FUN(J) .GE.FUN(I) )  GO  TO  16 
J  =  I 


16  CONTINUE 


c 

C  BASIC  LOOP.  ELIMINATE  EACH  WORST  VERTEX  IN  TURN. 

C  IT  MUST  BECOME  NO  LONGER  WORST,  NOT  MERELY  IMPROVED. 
C  FIND  NEXT -TO -WORST  VERTEX,  THE  ’ JN ’ TH  ONE. 

17  JN  =  1 

IF  (J.EQ. 1)  JN  =  2 
C 

DO  18  1  =  1,  K 
IF  (I.EQ.J)  GO  TO  18 
IF  (FUN(JN) .GE.FUN(I) )  GO  TO  18 
JN  =  I 

18  CONTINUE 
C 

C  LIMT  =  NUMBER  OF  MOVES  DURING  THIS  TRIAL  TOWARD  THE 
C  CENTROID  DUE  TO  FUNCTION  VALUE. 

LIMT  =  1 
C 

C  COMPUTE  CENTROID  AND  OVER  REFLECT  WORST  VERTEX. 

C 

DO  19  1=1, NV 
VT  =  V(I,J) 

SUM ( I )  =  SUM ( I ) - VT 

CEN(I)  =  SUM( I ) / FKM 

VT  =  BETA" CEN ( I ) -  ALPHA" VT 

IF  (IP.EQ.l)  VT  =  AINT (VT* . 5 ) 

C 

C  INSURE  THE  EXPLICIT  CONSTRAINTS  ARE  OBSERVED. 

19  V(I,J)  =  AMAX1 ( AMIN1 ( VT , BU ( I ) ) , BL ( I ) ) 

C 

NT  =  NT* 1 
C 

C  CHECK  FOR  IMPLICIT  CONSTRAINT  VIOLATION. 

C 

20  DO  25  N= 1 , NLIM 
NCE  =  NCE*  1 
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IF  (KE(V(L, J) ) .EQ.O)  GO  TO  26 


C 

C  EVERY  ' KV ‘ TH  TIME,  OVER-REFLECT  THE  OFFENDING  VERTEX 
C  THROUGH  THE  BEST  VERTEX. 

IF  (MOD ( N , KV ) . NE . 0  )  GO  TO  22 
CALL  FBV  ( K , FUN , M ) 

C 

DO  21  1=1, NV 

VT  =  BETA-'-V  ( I ,  M )  -  ALPHA "V  ( I ,  J  ) 

IF  (IP.EQ.l)  VT  =  AINT ( VT* . 5  ) 

21  V(I,J)  =  AMAX 1 ( AMIN  1 ( VT , BU ( I ) )  . 3L ( I  )  ) 

C 

GO  TO  24 
C 

C  CONSTRAINT  VIOLATION:  MOVE  NEW  POINT  TOWARD  CENTROID. 

C 

22  DO  23  1=1, NV 

VT  =  . 5  * ( CEN ( I ) *  V ( I ,  J  )  ) 

IF  (IP.EQ.l)  VT  =  AINT ( VT* . 5  ) 

V(I,J)  =  VT 

23  CONTINUE 
C 

24  NT  =  NT+ 1 

25  CONTINUE 
C 

IER  =  1 
C 

C  CANNOT  GET  FEASIBLE  VERTEX  BY  MOVING  TOWARD  CENTROID, 

C  OR  BY  OVER- REFLECTING  THRU  THE  BEST  VERTEX. 

IF  (NPR.LE.O)  GO  TO  42 
WRITE  (6,52)  NT , J 

CALL  BOUT  ( NT , NPT , NFE , NCE , NV , NVT , V . K , FUN , CEN , J ) 

GO  TO  42 
C 

C  FEASIBLE  VERTEX  FOUND,  EVALUATE  THE  OBJECTIVE  FUNCTION. 


26  NFE  =  NFE*  1 


FUNTRY  =  FE ( V ( 1 , J ) ) 

C 

C  TEST  TO  SEE  IF  FUNCTION  VALUE  HAS  NOT  CHANGED. 

AFO  =  ABS( FUNTRY -FUNOLD) 

AMX  =  AMAX 1  ( AB  S  (  E P  FUNOLD ),  E P  ) 

C 

C  ACTIVATE  THE  FOLLOWING  TWO  STATEMENTS 
C  FOR  DIAGNOSTIC  PURPOSES  ONLY. 

C  WRITE  (6,99)  J , AFO , AMX , FUNTRY , FUNOLD , FUN ( J ), FUN ( JN ) , 

C  INTFS , N 

C  99  FORMAT  ( IX, 13 , 6EI5 . 7 ,215 ) 

IF  (AFO. GT. AMX)  GO  TO  27 
NTFS  =  NTFS  +  I 
IF  (NTFS.LT.NCT)  GO  TO  28 
IER  =  0 

IF  (NPR.LE.O)  GO  TO  42 
WRITE  (6,53)  K 

CALL  BOUT  (NT , NPT , NFE , NCE , NV , NVT , V , K , FUN , CEN , 0 ) 

GO  TO  42 

27  NTFS  =  0 
C 

C  IS  THE  NEW  VERTEX  NO  LONGER  WORST? 

28  IF  (FUNTRY. LT. FUN ( JN ) )  GO  TO  34 
C 

C  TRIAL  VERTEX  IS  STILL  WORST;  ADJUST  TOWARD  CENTROID. 

C  EVERY  'KV'TH  TIME,  OVER-REFLECT  THE  OFFENDING  VERTEX 

C  THROUGH  THE  BEST  VERTEX. 

LIMT  =  LIMT+ I 

IF  (MOD (LIMT ,KV) . NE . 0 )  GO  TO  30 
CALL  FBV  ( K , FUN , M ) 

C 

DO  29  1=1, NV 

VT  =  BETA-’- V  ( I  ,  M )  -  ALPHA"  V  ( I  ,  J  ) 

IF  (IP.EQ.l)  VT  =  AINT ( VT +  .  5  ) 


108 


29  V(I,J)  =  AMAX 1 ( AMIN 1 ( VT , BU ( I ) ) , BL ( I ) ) 


C 

GO  TO  32 
C 

30  DO  31  1=1, NV 

VT  =  .5*(CEN(I)+V(I,J)) 

IF  (IP.EQ.l)  VT  =  AINT ( VT+ . 5 ) 

V(I,J)  =  VT 

31  CONTINUE 
C 

32  IF  (LIMT.LT.NLIM)  GO  TO  33 
C 

C  CANNOT  MAKE  THE  ' J ’ TH  VERTEX  NO  LONGER  WORST  BY 
C  DISPLACING  TOWARD  THE  CENTROID  OR  BY  OVER- REFLECTING 
C  THRU  THE  BEST  VERTEX. 

IER  =  2 

IF  (NPR  .LE.  0)  GO  TO  42 
WRITE  (6,52)  NT,  J 

CALL  BOUT  ( NT , NPT , NFE , NCE , NV , NVT , V , K , FUN , CEN , J ) 

GO  TO  42 

33  NT  =  NT+1  • 

GO  TO  20 

C 

C  SUCCESS:  WE  HAVE  A  REPLACEMENT  FOR  VERTEX  J. 

34  FUN ( J )  =  FUNTRY 
FUNOLD  =  FUNTRY 
NPT  =  NPT  1 

C 

C  STOP  AT  THE  100 'TH  PERMISSIBLE  TRIAL 
IF  (MOD (NPT , 100 ) . EQ . 0 )  GO  TO  48 
C 

DO  36  1=1, NV 
SUM ( I )  =  0. 

C 


109 


■  -V;  -> 


DO  35  N= 1 , K 


35  SUM ( I )  =  SUM (I)+V(I,N) 

C 

CEN(I)  =  SUM ( I ) / FK 

36  CONTINUE 
C 

LC  =  0 
GO  TO  39 
C 

37  DO  38  1=1, NV 

38  SUM (I )  =  SUM (I)+V(I,J) 

C 

LC  =  J 
C 

39  IF  (NPR.LE.O)  GO  TO  40 

IF  (MOD (NPT , NPR ) . NE . 0 )  GO  TO  40 
C 

CALL  BOUT  ( NT , NPT , NFE , NCE , NV , NVT , V , K , FUN , CEN , LC ) 

C 

C  HAS  THE  MAX.  NUMBER  OF  TRIALS  BEEN  REACHED  WITHOUT 
C  CONVERGENCE  ?  IF  NOT,  GO  TO  NEW  TRIAL. 

40  IF  (NT.GE.NTA)  GO  TO  41 
C 

C  NEXT -TO -WORST  VERTEX  NOW  BECOMES  WORST. 

J  =  JN 
GO  TO  17 

41  IER  =  3 

IF  (NPR.GT.O)  WRITE  (6,54) 

C 

C  COLLECTOR  POINT  FOR  ALL  ENDINGS. 


C  1)  CANNOT  DEVELOP  FEASIBLE  VERTEX.  IER  =  1 

C  2)  CANNOT  DEVELOP  A  NO-LONGER- WORST  VERTEX.  IER  =  2 

C  3)  FUNCTION  VALUE  UNCHANGED  FOR  K  TRIALS.  IER  =  0 

C  4)  LIMIT  ON  TRIALS  REACHED.  IER  =  3 

C  5)  CANNOT  FIND  FEASIBLE  VERTEX  AT  START.  IER  =-l 

42  CONTINUE 
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C  FIND  BEST  VERTEX. 

CALL  FBV  ( K , FUN , M ) 

IF  (IER.GE.3)  GO  TO  44 

C  RESTART  IF  THIS  SOLUTION  IS  SIGNIFICANTLY  BETTER  THAN 
C  THE  PREVIOUS  OR  IF  THIS  IS  THE  FIRST  TRY. 

IF  (NPR.LE.O)  GO  TO  43 
WRITE  (6,55)  (M, YMN , FUN (M) ) 

43  IF  (FUN(M) .GE. YMN)  GO  TO  47 

IF  (ABS(FUN(M)-YMN) . LE . AMAX1 (EP , EP*YMN) )  GO  TO  47 
C  GIVE  IT  ANOTHER  TRY  UNLESS  LIMIT  ON  TRIALS  REACHED. 

44  YMN  =  FUN (M) 

FUN ( 1 )  =  FUN (M ) 

DO  45  1=1, NV 
CEN(I)  =  V ( I , M) 

SUM(I)  =  V ( I , M ) 

45  V  ( 1 , 1 )  =  V  ( I ,  M ) 

DO  46  1=1, NVT 

46  XS(I)  =  V ( I , M ) 

IF  (IER.LT.3)  GO  TO  6 

47  IF  (NPR.LE.O)  GO  TO  48 

CALL  BOUT  (NT , NPT , NFE , NCE , NV , NVT , V , K , FUN ,V(1,M),-1) 
WRITE  (6,56)  FUN (M) 

48  RETURN 

49  FORMAT  ( 50H0INDEX  AND  DIRECTION  OF  OUTLYING 
1VARIABLE  AT  STARTI5 ) 

50  FORMAT  ( 50H0IMPLICIT  CONSTRAINT  VIOLATED  AT  START. 
1DEAD  END. ) 

51  FORMAT  (’OCANNOT  FIND  FEASIBLE ’, I 4 ,' TH  VERTEX  OR 
1CENTROID  AT  START.') 

52  FORMAT  ( 10H0AT  TRIAL  I4,54H  CANNOT  FIND  FEASIBLE 


1VERTEX  WHICH  IS  NO 

2LONGER  WORST, 14, 15X, 'RESTART  FROM  BEST  VERTEX.') 

53  FORMAT  (40H0 FUNCTION  HAS  BEEN  ALMOST  UNCHANGED 
IFOR  15, 7H  TRIALS) 

54  FORMAT  (27HOLIMIT  ON  TRIALS  EXCEEDED.  ) 

55  FORMAT  (’OBEST  VERTEX  IS  NO. ’,13,’  OLD  MIN  WAS 
1E15 . 7 ,  ’  NEW  MIN  IS  ’ ,E15 . 7) 

56  FORMAT  ( ’ OMIN  OBJECTIVE  FUNCTION  IS  ’,E15.7) 

END 

SUBROUTINE  FBV  (K , FUN ,M) 

DIMENSION  FUN ( 50 ) 

M  =  1 
C 

DO  1  1  =  2, K 

IF  (FUN(M) .LE.FUN(I) )  GO  TO  1 
M  =  I 
1  CONTINUE 
C 

RETURN 

END 

SUBROUTINE  BOUT  (NT , NPT , NFE , NCE , NV . NVT , V , K , FN , C , IK ) 
DIMENSION  V(50, 50) ,  FN(50),  C(25) 

WRITE  (6,4)  NT, NPT, NFE, NCE 
C 

DO  1  I  =  1 ,  K 

WRITE  (6,5)  FN ( I ) , (V(J,I),J=1, NV ) 

IF  (NVT.LE.NV)  GO  TO  1 
NVP  =  NV  *  1 

WRITE  (6,6)  (V(J,I) ,J=NVP,NVT) 

1  CONTINUE 
C 

IF  (IK.NE.O)  GO  TO  2 
C 

WRITE  (6,7)  (C(I) , I = 1 , NV ) 

RETURN 

LI2 


2  IF  (IK.GE.O)  GO  TO  3 
WRITE  (6,8)  (C(I),I=1,NV) 

RETURN 

3  WRITE  (6,9)  IK, (C(I) ,I=L,NV) 

RETURN 

4  FORMAT  ('ONO.  TOTAL  TRIALS  =  ’,I5,4X,'N0.  FEASIBLE 
ITRIALS  =  ’ , 

215 , 4X , ’ NO .  FUNCTION  EVALUATIONS  =  ' ,I5,4X,’N0. 

3 CONSTRAINT  EVALUATIONS 

4=  , 15 / ' 0  FUNCTION  VALUE’ ,6X, 'INDEPENDENT 

5VARIABLES/ DEPENDENT 

60R  IMPLICIT  CONSTRAINTS') 

5  FORMAT  (1H  , E18 . 7 , 2X , 7E14 . 7 / ( 2 IX , 7E14 . 7 ) ) 

6  FORMAT  (21X.7E14.7) 

7  FORMAT  (lOHOCENTROID  11X , 7E14 . 7/ (21X , 7E14 . 7 ) ) 

8  FORMAT  (’O  BEST  VERTEX ’ , 7X , 7E14 . 7 / ( 2 IX , 7E 14 . 7 ) ) 

9  FORMAT  ( ’ OCENTROID  LESS  VX ’ , 12 , 2X , 7E14 . 7 / ( 2 IX , 7E 14 . 7 ) ) 
END 

FUNCTION  FE (X ) 

DIMENSION  X ( 3 ) - 
COMMON  TDIFF 
CALL  PLANT (X) 

FE=TDIFF 

RETURN 

END 

FUNCTION  KE (X ) 

DIMENSION  X ( 3 ) 

KE  =  0 
RETURN 
'  END 
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